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The  increased  bitpot tance  of  software  for  embedded  avionics  systems  nas 
lea  to  an  increasing  desire  to  insure  tnat  avionics  software  meets 
very  strict  reliaoility  and  quality  goals,  however,  a  significant 
problem  in  assuring  suen  goals  are  met  is  tne  inability  of  Government 
personnel  to  accurately  predict  tne  reliability  of  an  avionics 
software  development  project.  Tnis  problem  has  been  expressed  at 
several  Government  and  industry  sponsored  conferences,  as  well  as  in 
documents  suen  as  the  Joint  Logistics  Commanders  Software  Reliability 
Working  Group  Report  (hovemoet  1975)  and  tne  Joint  Logistics 
Commanders  Software  Quality  Management  workshop  Report  (July  1979).  As 
a  result,  efforts  nave  been  initiated  to  develop  and  validate 
matnematical  models  for  predicting  tne  reliability  ana  error  content 
of  a  software  system,  however,  models  developed  to  date  have  not 
aaecuatley  aouresseu  tne  unique  features  of  avionics  software 
oevelopt  nents. 


This  effort  was  initiated  in  response  to  tne  need  for  developing 
software  reliability  prediction  mouels  applicable  to  avionics  software 
developments,  ana  fits  into  the  goals  of  RACC  IRC  ho.  5,  Software  Cost 
Reduction,  in  tne  subthrust  of  Software  Quality  (Software  Modeling) . 
This  report  summarizes  the  development  of  a  mathematical  model  for 
predicting  tne  reliability  ana  mean-time-to-failure  of  a  software 
development  unaer  the  assumptions  of  incomplete  information  available 
on  error  correction,  ana  discrete  versions  of  tne  software  being 
developed.  The  report  also  describes  tne  modified  nonlinear  search 
algorithm  developed  for  finding  model  parameters  ana  an  accompany ing 

vii 
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computer  program  for  operating  the  moael.  The  importance  of  this  model 
development  is  that  the  assumptions  underlying  this  model  more  closely 
reflect  the  actual  avionics  software  development  process  than  prior 
model  developments. 

The  theory  ana  model  algorithm  developed  under  this  effort  will  lead 
to  much  needed  predictive  measures  for  use  by  software  managers  of 
avionics  software  developments  in  adequately  tracking  those 
developnents  in  terms  of  reliability  and  mean-time-to-failure 
objectives.  More  importantly,  the  measures  developed  under  this  effort 
will  be  applicable  to  current  avionics  software  developments  and  thus 
help  to  produce  the  high  quality ,  low  cost  avionics  software  needed 
for  today's  aircraft. 

QJL.nlM' 

ALAN  N.  SUKERT 
Project  Engineer 
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1.0  Introduction 


1.1  Problem  Statement 

As  the  cost  and  complexity  of  computer  software  continue  to 
increase,  there  is  a  growing  need  for  accurate  determination  of  soft¬ 
ware  reliability.  Before  a  software  package  is  put  into  operation, 
there  is  a  testing  period  during  which  errors  are  detected  and  corrected. 
The  problem  with  which  we  are  concerned  is  the  estimation  of  certain 
reliability  parameters  from  the  error  data  generated  during  the  test 
phase.  Specifically,  we  wish  to  estimate  the  number  of  errors  remaining 
in  the  software  package  at  any  time,  and  the  mean  time  to  failure  (MTTF). 
Accurate  determination  of  these  parameters  could  reduce  the  cost  associ¬ 
ated  with  excessive  testing,  and  could  increase  the  confidence  with  which 
the  package  is  used. 

In  order  to  estimate  software  reliability,  it  is  necessary  to 
develop  an  appropriate  model  describing  the  error  detection  and  correc¬ 
tion  processes,  and  to  develop  procedur-es"For  estimating  the  parameters 
of  this  model  from  observed  error  data.  Our  intention  is  to  generalize 
certain  models  which  have  previously  been  used  for  this  purpose  in  order 
to  depict  more  accurately  an  actual  testing  environment.  In  addition, 
we  will  consider  a  somewhat  different  approach  to  the  estimation  of  the 
parameters  of  this  generalized  model. 

1.2  Previous  Work 

A  substantial  body  of  work  now  exists  on  the  application  of 
statistical  modeling  and  estimation  techniques  to  the  determination  of 
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software  reliability.  We  make  no  attempt  to  describe  all  this  work, 
but  rather  restrict  ourselves  to  those  efforts  which  are  directly  re¬ 
lated  to  our  own.  For  a  comprehensive  review  and  bibliography,  see  [1] 
or  [2]. 

One  of  the  most  widely-used  error  models  was  developed  by  Jelinski 
and  Moranda  [3].  A  similar  model  has  been  considered  by  Shooman  [4]  a 
others.  The  assumptions  about  the  error-detection  and  error-correction 
processes  which  underlie  this  model  are  the  following: 

(a)  The  error-detection  process  is  a  Poisson  process  whose 
detection  rate  is  constant  between  error  detections. 

(b)  The  error-detection  rate  at  the  time  prior  to  the  detection 
of  the  ith  error  is  a  function  of  i;  it  is  denoted  by  z±. 

It  is  commonly  assumed  that  is  proportional  to  the  number 
of  errors  in  the  program  at  detection  time.  This  can  be 
written  as: 

t  (Nq  -  i  +  i)  (1.1) 

where  Nq  is  the  initial  number  of  errors  and  <J>  is  a  positive 
constant.  An  alternative  assumption  is  that  the  detection 
rate  forms  a  geometric  progression 

=  Xa1  (1.2) 

with  both  X  and  a  being  positive  constants.  It  should  be 
noted  that  the  main  justification  for  (1.2)  is  the  improved 


convergence  of  Che  resulting  estimator  equations  [5], 

(c)  Error  detection  is  followed  by  an  immediate  correction. 
Consequently,  upon  detection  of  the  iC^  error,  the  number 
of  remaining  errors  drops  to  -  ij. 

(d)  The  debugging  process  is  perfect  and  no  new  errors  are 
generated  by  the  correction  process. 

These  assumptions,  although  restrictive,  were  initially  adopted 
by  most  of  the  researchers  in  the  field.  The  estimation  of  the  reli¬ 
ability  parameters  was  based  on  the  above  assumptions,  and  employed  the 
maximum  likelihood  (ML,  criterion  to  derive  the  best  estimates. 

It  is  realized  now  that  the  assumptions  given  above  are  quite 
restrictive  and  unrealistic  in  most  cases,  and  steps  have  been  taken  to 
make  the  model  more  realistic.  The  model  assumptions  have  been  changed 
to  comply  more  closely  with  the  real  process. 

Goel  [6]  has  considered  a  nonideal  debugging  process  in  which 
the  probability  of  correcting  an  error  is  p.  Based  on  this  assumption, 
an  analysis  of  the  resulting  model  is  performed.  Further  generalization 
is  suggested  by  Shooman  [7],  who  modified  both  assumptions  c  and  d  above 
concerning  the  error-correction  process.  According  to  the  modified 
model  of  Shooman,  the  correction  process  does  not  necessarily  proceed 
identically  to  the  detection  process,  and  new  errors  may  be  introduced. 
Denote  by  r^Ct),  rc(t),  and  r^(t)  the  rates  of  error  detection,  correc¬ 
tion,  and  new  error  generation,  respectively.  The  models  suggested  by 
Shooman  assume  different  relationships  between  these  rates.  The  main 
models  are: 
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Model  1 


rc(t)  -  6  rd(t)  (1.3) 

r  (t)  *  a  r  (t)  (1.4) 

D 

rc(t)  =  b  rd(t)  (1.5) 

rg(t)  =  a  n(t)  rd(t)  (1.6) 

where  n(t)  represents  the  number  of  errors  in  the  program. 

These  models  and  others  have  been  studied  by  Shooman,  and  the  re¬ 
sults  are  described  [7]. 

Another  generalization  of  the  original  model  concerns  the  assump¬ 
tion  that  the  corrections  are  implemented  continuously.  This  is  not 
consistent  with  actual  practice  in  which  a  program  is  replaced  by  a 
newer  version  at  discrete  times.  Between  the  replacement  times,  the 
program  undergoing  the  test  is  the  same  and  the  number  of  errors  in  it 
is  constant.  A  possible  solution  for  this  discrepancy  is  that  redis¬ 
covery  of  errors  should  not  be  counted.  However,  this  requires  the 
analysis  of  the  source  of  errors  in  order  to  determine  whether  the  er¬ 
ror  sources  are  the  same,  and  this  is  not  always  practical.  A  modified 
model  in  which  this  generalization  was  implemented  was  discussed  by 


and 

Model  2 
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Tal  [5]  and  by  Sukert  [2],  and  estimator  equations  for  use  with  this 
model  were  developed. 

These  generalizations,  along  with  some  additional  ones,  will  be 
incorporated  into  a  new  model.  The  new  model,  we  believe,  more  ac¬ 
curately  describes  an  actual  testing  environment.  We  will  first  dis¬ 
cuss  the  behavior  of  this  model  as  a  function  of  its  parameters  under 
the  simplifying  assumption  that  the  error  processes  are  deterministic 
rather  than  random. 

After  presenting  certain  results  for  deterministic  processes,  we 
will  then  show  results  using  simulated  random  error  data.  We  have 
developed  a  least-squares  search  procedure  for  estimating  the  model 
parameters,  and  will  discuss  its  convergence  behavior.  Recommendations 
are  made  toward  increased  utility,  and  toward  closer  coupling  of  the 
algorithm  to  information  in  real  test  data. 

The  true  test  of  the  usefulness  of  the  model  will  lie  in  its 
ability  to  describe  real  software  tests.  Thus,  there  remains  for  sub¬ 
sequent  work  the  application  of  the  model  to  enough  real  cases  to  draw 
conclusions  concerning  validity. 

One  of  the  difficulties  encountered  by  researchers  in  the  past 
has  been  the  inadequacy,  incompleteness,  and  ambiguity  of  available  test 
data.  We  found  some  of  these  same  problems  with  the  data  available  to 
us  during  this  work.  Hence,  we  include  coranents  regarding  data  require¬ 


ments. 


2.0  Development  and  Analysis  of  the  Model 


2.1  Assumptions 


In  order  to  develop  a  generalized  model  to  describe  the  error 
detection  and  correction  processes,  we  make  the  following  assumptions: 

(a)  The  error  detection  process  is  a  Poisson  process  whose  aver¬ 
age  rate  of  occurrence  is  proportional  at  any  time  to  the 
number  of  errors  present  in  the  software  package.  Denoting 
the  number  of  errors  present  at  time  t  by  N(t)  and  the  aver¬ 
age  error  occurrence  rate  by  rd(t),  we  have 


rd(t)  =  $N(t) 


where  4>  is  a  fixed  constant  of  proportionality. 

(b)  No  attempt  is  made  to  correct  detected  errors  at  the  time 
of  detection.  Instead,  a  new  and  corrected  version  of  the 
program  is  provided  to  the  testing  group  at  discrete  ("tape 
replacement")  times  t^,  t^,  • ....  Thus,  the  number 
of  errors  present  in  the  program  at  time  t,  t  <  t  <  t^+^, 

is  constant  and  equal  to  N^t^.  This  is  illustrated  in  Fig.  1. 

(c)  Of  the  detected  errors  reported  to  the  correcting  group,  some  are 
corrected  and  some  are  not.  In  addition,  new  errors  are 
generated.  Denote  the  cumulative  number  of  errors  corrected 

to  time  t  in  the  program  being  tested  by  Nc(t),  and  the 
cumulative  number  of  newly-generated  errors  by  Ng(t).  Both 
Nc  and  N  are  piecewise  constant  because  of  the  assumption 
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Fig.  1.  Number  of  errors  in  program  as  a  function  of  time. 
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that  a  new  version  of  the  program  is  provided  only  at  the 
discrete  times  t^,  t2>  ....  A  key  feature  of  our  model  is 
that  many  errors  may  be  detected,  corrected,  and  generated 
between  update  times.  At  any  time  t  during  testing,  we 
have 

N(t)  =  N  -  N  (t)  +  N  (t)  (2.2) 

o  c  g 

where  is  the  initial  number  of  errors  in  the  program. 

(d)  The  error  correction  rate  rc(t)  depends  on  both  the  error 
detection  rate  rd(t)  and  the  error  backlog  N^(t),  where 

Nb(t)  *  Nd(t)  -  Nc(t).  (2.3) 

For  simplicity,  we  assume  a  linear  relationship 

rc(t)  -  ard(t)  +  BN^t).  (2.4) 

The  addition  of  the  second  term  in  (2.4)  represents  a 
generalization  of  the  model  of  Shooman  [7]. 

(e)  The  rate  of  generation  of  new  errors  is  proportional  to  the 
error-correction  rate: 

rg(t)  -  yrc(t).  (2.5) 

(f)  The  error-detection  process  Nd(t)  is  precisely  known,  but 
the  error-correction  process  Nc(t)  is  unknown.  This  ap¬ 
pears  to  be  a  realistic  assumption  in  view  of  the  way 


error  correction  is  actually  performed.  The  error 
generation  process  is  also  unknown. 

In  the  above  we  tacitly  equate  software  "failure"  to  coding 
"faults".  In  effect,  we  include  in  a  and  $  the  proportionality  be¬ 
tween  the  two,  and  call  them  both  "errors". 

2.2  The  Model 

The  model  which  we  develop  is  actually  a  deterministic  model 
which  relates  the  expected  values  of  the  various  random  processes  in¬ 
volved.  The  required  connection  between  the  observed  sample  functions 
of  the  random  processes  involved  and  the  deterministic  model  is  estab¬ 
lished  by  means  of  an  estimation  algorithm  which  operates  on  the  ob¬ 
served  data  to  estimate  model  parameters.  The  deterministic  model 
will  be  described  first,  followed  by  a  discussion  of  the  estimation 
procedure. 

Taking  expected  values  of  (2.1)-(2.4)  yields  the  equations 


rd(t)  -  4>n(t) ,  (2.6) 

r£(t)  =■  ard(t)  +  Bn^Ct),  (2.7) 

n(t)  «  Nq  -  n£(t)  +  ng(t),  (2.8) 

n^t)  -  nd(t)  -  nc(t),  (2.9) 


ng(t)  -  Ync(t), 


(2.10) 


where  a  lower-case  n  denotes  the  expected  value  of  the  process  repre¬ 
sented  by  the  corresponding  upper-case  N. 

It  follows  from  the  relationship  between  r,  (t)  and  n,(t)  that 

u  d 


nd(t) 


t 

rd  (t)du. 

r\ 


(2.11) 


Similarly, 


ft 

nc(t)  =  I  rc(u)du.  (2.12) 

J  o 

The  model  represented  by  the  above  equations  can  be  viewed  as  a 
linear  system  with  sampling  and  feedback  as  shown  in  Fig.  2.  Our 
problem  is  now  one  of  system  identification:  Given  Nd(t),  estimate 
the  parameters  of  the  system  shown  in  Fig.  2.  Revisions  to  the  soft¬ 
ware  are  applied  at  time  instants  t^»  between  which  times  n(t)  remains 
constant.  The  system  therefore  is  treated  as  a  discrete-time  system. 
We  employ  the  usual  notation  k  in  place  of  the  argument  t^. 

The  four  system  equations  (2. 6-2. 9)  can  be  reduced  to  two: 


rd(k)  *  $  [Nq  -  (1  -  y)  nc(k)]  (2.13) 

rc(k)  -  ap  j^NQ  -  (1  -  y)  nc(k)j  +  6  [nd00  “  t»c(k)]  (2.14) 
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Fig.  2.  Block  diagram  representation  of  the  proposed 
model  for  the  error-detection  and  the  error- 
correction  processes. 

and  the  number  of  parameters  reduced  to  four: 

rd(k)  =  ♦a  [Na  ‘  nc(k)]  (2'15) 

rc(k)  -  ct$a  -  nc(k)J  +  6  |nd(k)  "  nc(k)]  (2.16) 

where 

♦  #  «  (1  -  yH,  \  -  No/(l  -  y).  (2.17) 

The  application  of  Laplace  transform  techniques  and  some  algebraic 
manipulation  similarly  lead  to  the  equivalent  block  diagram  shown  in 
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Fig.  3.  Simplified  model. 


Fig.  3.  The  identification  problem  reduces  to  the  estimation  of  the 
four  parameters  N  ,  4>  >  a,  and  6.  Note  that  N  is  the  sum  of  the 

3  3  3 

initial  errors  N  and  all  errors  which  are  subsequently  generated 
o 

during  the  correction  process.  Note  further  that  the  ultimately 


(2.19) 


nd(k) 


n  (k) 
k  c 


and  using  the  usual  approximation,  which  in  our  case  is  exact, 


r(k)  *  -  n 'k)  •  T(k)  -  t(kfl)  -  t(k) 


(2.20) 


the  model  becomes 


n(k+l)  *  Ln(k)  +  N  <f>  B 
a  a 


(2.21) 


where 


-*T(k) 

a 


ST (k)  1  -  T(k)[s  +  a<t.a]y  \cxT(k) 


(2.22) 


When  tape  replacement  occurs  at  uniform  time  intervals,  T  is  constant 
over  k  and  the  system  is  seen  to  be  stationary,  and  the  equations  can  be 
solved  immediately  by  successive  evaluation: 


n(l)  -  H  ♦  B,  n(0)  -  0 

a  a 

n(2)  -  N  *  (L  +  1)B 

a  a 

n(3)  -  M  ♦  (l2  L  +  l)  B 


-  k_1  1 

n(k)  -  H  ♦  l  L-'B 

a  4  j-0 
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and  applying  the  familiar  procedure  for  the  geometric  sum, 

Ln(k)  -  n(k)  *  N  ♦  ( L^B  -  b) 
a  a 

which  gives  for  the  state  at  the  k  tape  replacement  time, 

n(k)  =  Na<fra  (L  -  I)-1  (,Lk  -  i)  B  (2.23) 

The  increment  ^(k)  =  n(k)  -  n(k-l)  at  the  kC^  tape  replacement  time  is 
given  by: 

!(k)  =  NJ  (L  -  1)_1  (Lk  -  Lk_1)  B 

a  8i 

*  NJ  Lk_1  B  (2.24) 

a  a 

2. 3  Model  Behavior 

Note  from  the  discrete  state  equations  above  that  the  parameter 
Na  is  simply  a  scale  factor  on  the  state  n.  Recall  also  that  the 
initial  slope  of  n ,  (k)  is  N  $  ,  and  that  of  n  (k)  is  aN  <j>  ,  regardless 
of  the  value  of  0.  Furthermore,  for  6*0,  n^(k)  and  n£(k)  maintain 
the  constant  ratio  nc(k)/n{j(k)  ■  a  <  1,  and,  of  course,  coincide  as 
a  ■*  1. 

The  effect  of  6  >  0  is  to  increase  the  error  correction  rate, 
and  therefore  increase  nc(k),  especially  for  the  larger  differences 
n^(k)  -  i»c(k)  (backlog)  which  tend  to  occur  later  in  the  test  program. 

The  resulting  decrease  in  remaining  errors  Nfl  -  n£(k)  causes  the  detected 
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error  curve  n^k)  to  be  bent  downward.  Thus  the  effect  of  8  is  to  draw 
the  two  curves  together.  Figures  4  and  5  display  this  effect  for  0  <  8 
<  0.5.  The  "bending"  of  the  curves  due  to  8,  together  with  the  effects 
of  the  discrete  nature  of  the  model,  are  expected  to  occur  in  real  data. 

2.4  The  Data  Simulator 

RELY  I  contains  a  data  simulator  for  the  purposes  of  study  and 
experimentation.  The  simulator  is  an  optional  source  of  input  data  to 
the  estimator  (see  Appendix  C).  The  simulator  reads  from  input  cards 
the  nominal  parameter  values,  a,  8,  <p  ,  N  ,  the  time  interval  T,  the 

cL 

number  of  test  intervals  K,  and  an  input  initial  random  number  (RRR) , 

and  computes  the  associated  software  test  history  A  (k).  The  random 

number  RRR  is  changed  by  the  investigator  when  he  wishes  a  different 

sample  of  the  random  data  set  A, (k)  (see  Appendix  A,  RNDTA,  RANDEX) . 

a 
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0-0.5 


Cumulative  detected  errors  n ,  (k)  for  0  <  8  <  0.5 


I 


3.0  The  Estimation  Algorithm 

3.1  The  Estimation  Problem  and  Method 

Having  described  the  model,  we  turn  to  the  parameter  estimation 
algorithm  which  estimates  the  values  of  the  model  parameters  correspond¬ 
ing  to  a  given  set  of  real  test  data.  The  resulting  parameter  estimates 
provide  the  reliability  information  sought  regarding  the  tested  software 
package. 

Though  the  model  is  linear  in  the  sense  that  the  equations  are 

linear  in  the  state  n,  the  model  equations  are  nevertheless  nonlinear 

in  the  parameters  £  (i.e.,  in  a,  8,  4>  ,  N  ).  Determining  the  parameter 

3  3i 

— ► 

values  corresponding  to  a  given  set  of  real  test  data  A(k)  is  then 
a  nonlinear  estimation  problem. 

Nonlinear  parameter  estimation  methods,  in  general,  are  iterative 
procedures  in  which  the  estimate  is  approached  from  some  initial  guess 
for  the  parameter  values,  in  steps  which  successively  decrease  a  cost 
functional  J.  Since  our  purpose  is  to  determine  the  parameter  values  0 
for  which  the  solution  {^(k)  of  the  model  equations  approximates 
the  measured  function  (or  sequence)  A^(k),  we  choose  the  cost  functional 


J  to  be  the  sum  of  the  squares  of  the  residuals,  6,00  -  A , (k) ,  viz., 

d  a 


[5d(k)  -  Vk>] 


(3.1) 


Minimizing  this  cost  functional,  then,  minimizes  the  difference  between 


the  observed  function  A,(k)  and  its  expected  value  6,(k)  in  the  least 

a  a 


squares  sense. 


-  3.1  - 


\  >  9  '»»» 


Of  the  numerous  methods  described  in  the  literature,  both  direct 
search  methods  (Fletcher  [9])  and  gradient  methods  (Bard  [8]),  the 
gradient  methods  are  generally  preferred  when  the  computation  of  the 
derivatives  of  J  is  not  prohibitive.  Gradient  methods,  in  principle, 
step  from  one  point  6  in  parameter  space  to  the  next  according 

to 


TiRi8i 


(3.2) 


where  is  the  gradient  of  J  evaluated  at  0^,  R^  is  some  matrix  which 

til  “► 

operates  on  the  gradient  to  define  the  i  step  direction  R^g^,  and 
is  a  scalar  which  determines  the  step  size.  The  methods  differ  in  what 
each  employs  for  R^,  i.e.,  in  the  step  direction  each  takes  relative  to 
the  gradient.  The  method  of  steepest  descent,  for  example,  uses  the 
identity  matrix  for  R^,  so  that  the  step  direction  is  opposite  to  that 
of  the  gradient.  This  is  "the  steepest  way  down"  locally  but  tends  to  be 
less  efficient  and  therefore  less  desirable  than  methods  which  use  second 
order  information  about  the  surface  J(8). 

The  Newton-Raphson  method  uses  for  R^  the  inverse  of  the  Hessian, 
the  matrix  of  the  second  partial  derivatives, 


H 


mn 


a2j 


36  86  ’ 
m  n 


(3.3) 


of  the  cost  functional. 

i 

i 


Notice  that  the  Taylor  series  expansion  of  J  to  second  order  terms. 
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J  '  J1  *  P  -  5t)  +  7  (®  '  «i)T  Hi  (5 


has  an  extremum. 


-  -  *1 + h  ?  -  h)  -  “■ 


0  *  0^  -  (H^  nonsingular) 


so  if  Ri  *  anc*  is  Quadratic,  then  ei+1  coincides  with  the 

extremum.  The  Newton-Raphson  method  in  this  case  converges  in  a  single 
iteration.  This  method  is  quite  efficient  even  for  nonquadratic  J,  but 
only  if  H  is  positive  definite.  This  latter  condition  is  the  principal 
weakness  of  the  method.  The  Marquardt  method  meets  this  weakness  by  guarantee¬ 
ing  positive  definiteness  in  R^  by  adding  to  (or  to  some  convenient 

2 

approximation  of  H^)  a  variable  amount  of  a  positive  definite  matrix  C^: 


*i  ■  (Hi +  vi) 


and  suggests  be  a  matrix  of  the  diagonal  elements  of  H^,  viz.. 


C.  -  H. 
i,ss  1  i,ss' 


For  sufficiently  large  then  is  positive  definite,  even  when  H1 

is  not.  The  Marquardt  method  behaves  as  the  Newton-Raphson  for  small 
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A^,  but  where  larger  A^  is  necessary  it  steps  nevertheless  in  some  accept 
able  (downward)  direction,  A  step  is  said  to  be  acceptable  if  it  de¬ 
creases  J.  If  A  is  large  and  has  low  condition  number  (eigenvalues 
of  near-equal  magnitude),  the  method  approximates  that  of  steepest 
descent.  The  Marquardt  method  varies  from  step  to  step  according  to  A^, 
between  the  behavior  of  the  Newton-Raphson  method  and  that  of  steepest 
descent . 


3.2  Description  of  the  Search  Algorithm 

The  program,  RELY  I,  uses  the  above  Marquardt  R^,  i.e., 

Vl  ’  91  -  'i(Hl  +  Vi)  *i  <3'6> 

and  selects  r  or  A^  from  step  to  step  according  to  the  procedure 
described  below.  Essentially  the  program  progresses  in  one  or  the  other  of 
two  modes.  In  mode  A,  A^  is  fixed  while  the  largest  (0.0001  <  <  1) 

is  sought  which  results  in  an  acceptable  step  size.  If  the  sought  is 
found,  the  program  continues  in  mode  A  preferring  smaller  and  smaller 
values  of  A  (more  nearly  Newton-Raphson).  If  at  any  point  insufficient 
progress  is  being  made  in  mode  A,  the  routine  moves  to  mode  B,  in  which 
is  initially  fixed,  and  A  is  successively  increased  until  an  accept¬ 
able  step  direction  is  reached.  In  mode  B,  when  a  sufficiently  large  A 
is  reached,  then  the  program  steps  In  that  direction  until  J  begins  to 
increase,  or  until,  for  large  J,  J  has  decreased  more  than  10  percent, 
at  which  point  the  routine  returns  to  mode  A.  In  short,  when  progress 
is  slow  in  mode  A,  the  program  resorts  to  mode  B  to  move  to  a  different 
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"locality".  "Progress"  in  mode  B  is  deliberately  restricted  for  large 
J  due  to  experience  which  indicates  that  mode  B  for  large  J  tends  to 
settle  into  local  minima.  The  program  terminates  when  J  becomes  less 
than  a  predetermined  value  (ERR),  or  upon  a  time  limit  for  machine 
computation. 

More  specifically,  the  estimator  proceeds  as  follows:  Given 
initial  guess  0  and  A^  =  1,  i  «=  0: 

Mode  A 

1.  Compute  cost  J.  and  step  direction  R^g^. 

2.  If  <  ERR  terminate,  otherwise  determine  an  acceptable 

step  size  in  the  following  way: 

a.  Compute  a  x^  such  that  twice  the  associated  step 
causes  6  =  5;  i.e., 

Ti-  (5-  - 

b.  If  such  a  step  causes  <J>  >0.2,  choose  instead 

a 

Ti  -  (°-2  -  »al)/[2l’Wl 

c.  If  the  resulting  x^  >  1,  set  x^  =  1. 

d.  If  the  resulting  x^  <  .0001,  jump  to  mode  B. 

e.  Limit  0  <  3  <  10.  Compute  J^+^. 

f.  If  ^i+1  -  jumP  t0  item  4  below. 

3.  Accept  and  reduce  A;  i.e., 

a.  Set  0i  -  61+1,  At  -  A^/10. 
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b.  If  J  has  decreased  less  than  1  percent  in  more  than 
five  iterations  (reductions  of  X  in  item  3. a)  since 
passing  through  mode  B,  jump  to  mode  B.  Otherwise 
continue  in  mode  A  (jump  to  item  1  above). 

4.  Reduce  x^  by  a  factor  of  10.  If  the  resulting  J^+1  <  jump 

to  item  3  above,  otherwise  repeat  item  4  above  up  to  five 
times  (according  to  counter  INDEX) .  If  J  does  not  decrease 
with  five  reductions  of  x^,  Jump  to  mode  B.  If  <  ERR, 

terminate. 

Mode  B 

5.  Fix  x  =  0.1,  set  X  »  0.01. 

6.  Increase  X  by  a  factor  of  10,  increment  the  count  I CLAM, 
determine  the  corresponding  step  direction  +  X^C^ 


parameter  set  and  Ji+^.  If  Ji+1  -  rePeat  item  6. 


7.  Accept  0i+1  (i.e.,  set  9i 


9i+l^  anc*  set  JX 


l+l’ 


8.  Increase  x  by  a  factor  of  5  ,  i 
1  S-l 


I CLAM. 


9*  Try  Vi =  ®i  -  Ti(Hi +  x±cl )  «i*  if  Ji+i ”  Ji or  if  Ji  ” 

50  and  J  ./.L  <  0.9,  accept  3.  and  return  to  item  1  above, 

lri  A  X 

otherwise  accept  -*■  8^  anc*  rePeat  ftem  8. 

The  algorithm  is  depicted  in  the  flow  diagram  of  Fig.  6. 
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3.7 


Experience  during  development  of  RELY  I  proved  a  to  be  insuffi¬ 
ciently  independent  of  the  other  parameters  to  warrant  a  fourth  degree 
of  freedom  in  the  search  process.  The  computer  program  therefore  was 
modified  to  accept  a  priori  the  estimate  of  a,  and  to  search  in  three 
dimensions  for  the  values  of  6,  <j>  ,  and  N  .  From  these  the  estimated 

H  3. 

number  of  remaining  errors, 

na  =  na(k)  =  Na  "  nc&>»  k  =  1.  2>  3,  •••  K 
and  mean  time  to  failure 

MTTF  = 

n 

a  a 

are  computed.  The  latter  two  computed  quantities,  the  sought  software 
reliability  factors,  were  found  to  be  essentially  insensitive  to  reason¬ 
able  a  priori  estimates  of  a.  Results  below  include  cases  of  correct 
and  incorrect  fixed  a. 

3.3  Estimation  Results 

Results  are  tabulated  and  displayed  in  histograms  below  for 
several  simulated  random  data  examples.  Examples  I  and  II  differ  in 
the  selection  of  K  and  T  to  vary  the  number  of  errors,  N  ,  remaining  in 
the  software.  Example  I  uses  test  interval  length  T  ■  1.5  and  60  inter¬ 
vals  which  leaves  about  30  remaining  of  the  initial  200  software  errors. 
Example  II  uses  longer  test  Intervals,  T  -  8,  and  fewer  intervals,  K  - 
20,  to  leave  about  5  errors  remaining.  Example  III  corresponds  to  a 


larger  software  system  having  a  considerably  larger  number  of  initial 
errors  (N  *  1000),  smaller  error  detection  rate  (<J>  =  0.01),  but  a 

d  3. 

better  correction  rate  (a  =  0.8),  and  the  longer  test  intervals  (T  * 

8.0).  Finally,  the  fourth  example  demonstrates  the  insensitivity  to  the 

fixed  value  of  a.  Example  IV  essentially  is  Example  I  with  a  fixed  at 

0.5  instead  of  the  "true"  value,  0.7. 

Before  examining  the  results,  we  anticipate  the  nature  of  the 

distributions  by  analytically  determining  the  mean  and  variance  for 

the  simple  single  interval  (K  =  1)  case.  Let  the  observed  number  of 

detected  errors  N,  be  Poisson  with  mean  and  variance  pN  ,  where  p  =*  <t>  T . 

d  a  a 

We  minimize  the  squared  error  J, 


J  -  (»»,  -  Nd)2 

M-  .  2p  (»H.  -  »„)  -  0 
a  v 


to  obtain  an  estimate 


N 


a 


P 


which  has  mean 


and  variance 
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2  2  2  2 
pN  +  p  S  -  p  N 

a _ a _ a 

2 

P 


N 

a 

P 


The  number  of  remaining  errors. 


is  estimated 


N  =  N  -  N  -  N  P-") 

R  a  d  d  ^  p  / 

with  unbiased  mean 


E  !S)  •  E  {Sa  -  Nd}  -  -  *d 

and  with  root  mean  squared  difference  from  its  true  value. 
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Notice  that  the  value  of  the  latter  quantity  corresponding  to: 

a.  Example  I:  Let  N0  =  N  -  pN  ,  or  (p  *  1  -  31/200  =  .845)  is 

K  a  a 


200 

.845 


15 


b.  Example  II:  (p  ■  1  -  5/200  -  0.975)  is 


200 


H  .975 


14 


c.  Example  III:  (p  -  1  -  185/1000  =  0.815)  is 


1000 

.815 


35 


One  would  expect  these  values  to  approximate  the  standard  deviations 
a  (Nr)  for  the  respective  multiple-interval  cases  (though  perhaps  with 
less  validity  when  NR/Nfl  is  small).  The  cj(Nr)  indicated  below  for  the 
four  examples  then  are  of  the  magnitude  to  be  expected.  Table  1  lists 
numerical  information  from  the  four  examples.  Examination  of  results 
from  the  four  simulated  examples  indicates  that  the  estimator  produces 
reasonable  estimates  of  the  reliability  parameters  NR  and  MTTF. 
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Table  1.  Estimation  results 


Item 

Example :  I 

II 

III 

IV 

"True"  Values 

a 

.700 

.700 

.800 

.700 

S 

.100 

.100 

.100 

.100 

*a 

.020 

.020 

.010 

.020 

N 

a 

200 

200 

1000 

200 

nr 

31.0 

4.93 

185 

30.6 

MTTF 

1.61 

10.1 

0.538 

1.64 

Time  Interval 

T 

1.5 

8.0 

8.0 

1.5 

Number  of  Intervals 

K 

60. 

20.0 

20.0 

60.0 

A  Priori  a 

af 

.700 

.700 

.800 

.500 

Initial  Guess 

6 

.300 

.300 

0 

0 

♦. 

.010 

.010 

.02 

.010 

N 

a 

300 

300 

1500 

300 

Estimated  Values 

h 

26.3 

5.03 

195 

29.7 

°  <V 

11.2 

3.07 

40.4 

10.8 

MTTF 

2.19 

12.6 

.543 

1.75 

a  (MTTF) 

.65 

7.43 

.093 

.33 
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Fig.  8 
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.  Histograms  of  reliability  parameters  for  Example  II. 
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4.0  Conclusions  and  Recommendations 


4.1  RELY  I 

We  have  developed  and  displayed  a  model  which  we  believe  more 
accurately  describes  an  actual  testing  environment  of  a  large  software 
package.  This  new  generalized  model  has  been  incorporated  in  an  estima¬ 
tion  algorithm  for  the  purpose  of  discerning  reliability  of  the  software 
from  its  test  data.  The  first  version  of  the  algorithm  RELY  I  described 
here  converges  in  a  given  region  of  interest  of  the  model  parameters. 

RELY  I  is  applicable  to  software  test  cases  where  tape  replacement  (soft¬ 
ware  revision)  occurs  at  uniform  intervals  of  time,  and  where  sufficiently 
reliable  information  is  available  concerning  the  number  of  errors  detected 
during  each  of  the  successive  intervals. 


Data  Requirements 


Data  required  for  RELY  I  are  simple,  viz.,  the  time  interval  T 

between  software  revisions  (tape  replacements),  and  the  sequence  A ,  (k) 

d 

during  each  of  the  K  successive  versions  of  the  software,  where  k  =  1, 

2,  ...,  K.  Secondly,  the  data  must  be  from  a  process  of  the  type  upon 
which  the  assumptions  of  the  model  were  based,  viz.,  the  testing  of 
large-scale  software  packages  such  as  that  in  the  F-16  control  system. 

There  must  be  an  identifiable  single  continuous  line  of  soft¬ 
ware  package  Identity  throughout  the  test  process.  The  package  passes 
successively  through  a  sequence  of  versions  k,  k  *  1,  2,  . . . ,  K.  At 
any  given  time  during  test,  the  software  is  in  only  one  of  the  versions. 
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1 


Each  version,  k,  is 


"all"  of  which  software  version  is  being  tested, 
identical  to  the  preceding  version,  k  -  1,  and  succeeding  version,  k  +  1, 
except  for  the  software  corrections  "counted"  in  Ac(k)  and  A£(k  +  1), 
respectively.  Figure  11  indicates  the  time  relationship  of  the  several  se¬ 
quential  quantities.  The  requirement  is  that  A^(k)  be  precisely  known 
for  each  version  k,  where  all  versions  are  identified  and  satisfy  this 
single  and  continuous  identity  as  described.  This  requirement  is 
violated  if  a  major  untested  version  is  suddenly  introduced  midstream, 
or  if  an  alternate  part  of  the  software  package  simultaneously  being 
tested  suddenly  is  adopted.  The  generated  error  feature  can  accommodate 
a  minor  amount  of  this  kind  of  violation,  but  generated  errors  are 
modeled  as  occurring  as  a  constant  proportion  of  the  correction  rate. 

Errors  are  usually  classified  into  certain  arbitrary  categories, 
ranging  from  those  obviously  to  be  counted,  to  those  of  doubtful 
pertinence  (obviously  "repeated"  errors,  errors  associated  purely  with 
erroneous  test  conduct,  etc.).  Suffice  it  here  to  suggest  that  the 
criterion  for  counting  a  given  error  or  not  will  be  related  to  its 
likelihood  of  occurrence,  and  its  interpretation  as  a  "failure",  under 
operational  conditions. 

A. 3  Recommendations 

Recommendations  toward  improved  interfacing  with  information  in  a 
real  test  process  (thus  taking  greater  advantage  of  inherent  features 

1  That  is,  all  the  parts  of  the  software  package  are  being  exercised  in 
a  manner  representing  that  for  which  the  reliability  factors,  e.g., 
MTTF,  are  to  be  applied  later. 
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Errors,  cumulative:  0 


Errors,  incremental:  0 


Version: 

Error  rate: 

Test  time: 
k: 


. 

(initial) 

J 

(1) 

1 

(2) 

r(0) 

^(1) 

r(2) 

Fig.  11.  Tested  software  in  only  one  version  at  a  time, 

identical  to  neighboring  versions  except  for  the 
changes  counted  in  A  at  the  respective  boundaries 


of  the  underlying  new  model)  include  the  following  further  work: 

1.  Revise  the  KOST  subroutine  to  accommodate  nonuniform  test 
intervals  T(k),  by  using  a  finite  difference  technique  for 
the  solution  of  ^(k,  ^).  The  increased  utility  is  expected 
to  far  outweigh  the  lesser  analytic  tractability  of  the 
resulting  system  and  the  possible  increase  in  required 
computation  time. 

2.  Apply  the  algorithm  to  real  data.  Available  data  should  be 
gathered,  studied,  and  adapted,  by  interpretation  and 
transformations,  to  the  requirements  of  RELY.  Residual 
functions  over  a  variety  of  cases  will  indicate  how  well 
the  model  represents  the  real  test  process.  Experience 
will  lead  to  further  recommendations  concerning  data  re¬ 
quirements,  and  to-  possible  improvements  in  RELY  such  as 
provisions  for  using  information  in  real  data  concerning 
error  correction  and  error  generation. 

For  example,  the  quantity  na(k). 


n  (k)  «  N  -  n  (k) 
a  a  c 


N  -  n  (k)  +  n  (k) 
.  °  c  S 

1  -  Y 


n(k) 
1  -  Y 


is  the  augmented  number  of  errors  remaining  in  the  tested 
software.  That  is,  nfl(k)  is  the  number  of  errors  which  would 
be  detected  henceforth  if  the  testing  process  were  to 
continue  indefinitely,  including  those  generated  after  time 
k.  The  number  of  errors  remaining  in  the  software,  exluding 
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those  yet  to  be  generated  in  the  correction  process,  is 


n(k)  =  (1  -  y)  n  (k) 
a 


The  parameter  y  is  assumed  not  observable  in  the  present 

implementation  of  the  model.  If,  however,  among  the 

detected  errors,  generated  errors  are  distinguishable  from 

original  errors,  then  the  additional  quantity  n  , (k) ,  the 

g“ 

number  of  generated  errors  detected,  is  available.  The  model 

state  is  easily  augmented  to  include  n  , (k) .  The  model  re- 

gd 

mains  unchanged  but,  to  the  extent  that  the  additional 
information  is  available  in  test  data,  the  model  parameter 
Y  becomes  observable. 

Experience  in  the  development  of  RELY  I  suggests  further 
investigation  of  the  nature  of  the  J(8)  surface.  Such  investigation 
should  include  also  the  surface  associated  with  the  alternative  cost 
functional  using  cumulative  functions  N(k),  rather  than  the  incremental 
number  of  errors  A(k).  Convergence  properties  in  certain  regions  of 
parameter  space  may  be  significantly  improved  using  N(k)  rather  than 
their  derivatives  A(k).  Indeed,  parallel  computation  using  each, 
respectively,  may  prove  both  feasible  and  advantageous.  Another 
gradient  type  parameter  estimation  method,  such  as  the  Fletcher-Powell 
deflected  gradient  method,  may  also  prove  more  efficient  with  the 
alternative  functional. 
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APPENDIX  A 


MAIN  AND  SUBROUTINE  DESCRIPTIONS 

1.  MAIN  and  Subroutine  Diagram 

The  subroutines  of  RELY  I  are  indicated  in  Fig.  A.l.  Internal 


Fig.  A.l.  RELY  I  subroutine  diagram. 

subroutines  EIGMIN  and  MARQ  are  shown,  as  well  as  the  UNIVAC  MATH-PACK 
library  subroutines  RANDEX,  TRIDMX,  EIGVAL,  and  EIGVEC.  A  brief 
description  of  MAIN  and  its  subroutines  follows. 


2.  MAIN 


MAIN  reads  input  data,  executes  the  estimation  algorithm  (see 

Sec.  3.2),  and  prints  output.  It  also  computes  simulated  random  data 

A,(k)  under  the  ISIM  *  1  option  (see  Appendix  C) . 
d 

i 


t 

! 


f 


3.  Subroutine  RNDTA  (A,  B,  P,  NA.  T.  V.  RR,  JJJ) 


From  a,  8,  ,  N  ,  T,  the  input  initial  random  numbers  and  L 

(corresponding  to  RNDTA  variables  A,  B,  P,  NA,  T,  RR,  JJJ,  respectively, 
and  corresponding  to  the  MAIN  variables  ALPH,  BETA,  PHI,  NA,  TD,  RRR, 

NN,  respectively),  RNDTA  computes  the  sequence  A, (k) ,  k  *  1,  2,  ...  L 

Q 

(the  RNDTA  variable  V,  and  MAIN  variable  S).  The  resulting  sequence 
A^(k)  is  used  as  simulated  data,  the  (incremental)  number  of  errors 
detected  successively  in  each  software  test  interval.  The  initial 
random  number  RRR  is  passed  through  POISS  to  RANDEX. 

RNDTA,  at  each  interval  k.  Integrates  the  system  equations: 

nd(k)  =  Nj(k  -  1)  +  <p&T  ^Na  -  nc  (k  -  1)) 

nc(k)  -  nc(k  -  1)  +  ot$aT  (Na  -  nc  (k  -  1))  +  8T[Nd(k  -  1)  -  n,(k  -  1)] 


Nd(0)  -  nc(0)  -  0 


using  the  cumulative  (random)  number  of  errors  detected  Nd(k  -  1),  to 
obtain  nc(k)  for  use  in  POISS.  Subroutine  POISS  generates  the  random 
integer  Ad(k)  according  to  the  mean  detection  rate  rd(k)  **  $fNa  “  nc 
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From  na  (*Na  -  nc),  $a»  T,  the  input  initial  random  number  RRR, 
and  k (POISS  variables  DD,  PP,  TT,  RRRR,  and  KKK,  respectively,  cor¬ 
responding  to  RNDTA  variables  D,  RP,  RT,  RQ,  and  K) ,  POISS  computes 
the  value  A^Ck),  according  to 

m 

l  C(i)  <  T,  A  (k)  -  m 
i-1  a 

The  random  sequence  C(i),  i  =  1,  2,  ...  100,  with  exponential  distribu- 
rdC 

tion  function  1  -  e  ,  r^  -  ^(Na  “  nc) ’  is  generated  by  RANDEX.  The 
starting  random  number  required  by  RANDEX  in  C(l)  is  the  input  initial 
random  number  RRR  for  k  *  1,  and  is  the  preceding  random  numbeT  r(i22'9 
for  k  >  1,  where  R  is  the  value  C(100)  previously  computed  for  the 


5.  Subroutine  RANDEX  (C,  100,  U) 

Reference:  UNIVAC  Large  Scale  System  MATH-PACK,  Programmer's 
Reference,  UP-7542,  Rev.  1. 

RANDEX  produces  a  set  of  100  pseudo-random  numbers  C  with  exponen¬ 
tial  distribution 


1 


C 


-UC 


by  operating  on  a  uniformly  distributed  variate  X,  according  to  the 
inverse  transform  method 


C 


-ln(l  -  X) 
U 


RANDEX  uses  two  other  UNIVAC  MATH-PACK  subroutines  RANDU  and  RANDN.  RANDU 

generates  X,  0  <  X  <  1,  for  which  computation  it  calls  RANDN  for  random 

35  35 

integers  0  <  I  <  2  .  RANDEX  requires  an  initial  value,  0  <  C(l)  <  2  , 

different  integer  parts  of  which  produce  different  output  sequences. 
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6.  Subroutine  KOST  (AA.  BB,  PP,  NNA,  SS,  NJJ.  TT. 

NMN ,  DT.  H.  GJ.  DND.  COST,  PNC,  RMTTF,  ZNC,  RERR) 

Given  values  of  the  parameters  time  interval  T,  number  of 
intervals  K,  test  data  A^Ck),  KOST  computes  the  incremental  error 
sequences  (see  Sec.  2.2) 


6  (k)  «  N  $  (1  0)L  B 

d  a  a 


6  (k)  -  N  <t>  (0  l)Lk-1  B 
c  ara 

where  (1  0)  and  (0  1)  are  the  transposes  of  the  vectors  and 

respectively,  and 


l'1  -*aT 


1  -  T(6  + 


KOST  further  computes  the  cost  scalar  (see  Sec.  3.1) 


j  -  l  kw  - * d (k)] 

k«=l  L 


the  gradient  vector  components 


..  K  3«,<k) 

W’W  Vk>  -if— 

m  k»l  m 


and  the  Hessian  matrix  elements 


a2  K  <k)  36.  (k)  96.  (k) 

90  96  .L  dVK'  90  90  90  90 

n  m  k“l  n  m  n  m 
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LjA. 


KOST  also  computes  the  associated  estimate  of  total  errors  corrected 


n 

c 


K 


l 

k=l 


5  00 

c 


the  number  of  errors  remaining 


n 


R 


N 

a 


n 

c 


and  the  mean  time  to  failure 


MTTF 


<p 

a  R 


1 

The  first  derivatives  above  are  given  by: 


5&d  00 
“se 

m 


N  <J>  (1 
a  a 


0)L 


k-2 


(k-D  B  4 

m 


L 


3B 

36 

m 


+  6d<k> 


,  3<J>  ,  3N 

1  a  _1 _ a 

$  39  +  N  30 

am  am 


where 


36i  fl 

i  *  m 

39m  jo 

i  ^  m 

and 
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The  second  derivatives  are: 


3  6d(k) 

30  36 

n  m 


N  <J>  (1 
a  a 


i  (k-2)Lk~3  ~  |(k-l)  ~  B  +  L 


16,00 

I  d 


1 

3*a 

i 

3Na 

<t> 

30 

N 

36 

a 

m 

a 

m 

-  A. 8 


1 .  Subroutine  EIGMIN  (HH,  CORK,  GGJ.  DDAJ1,  EH,  EV) 

Given  the  Hessian  HH(4,  4)  and  the  gradient  GGJ(4)  of  the  cost 
functional  J(6),  and  the  current  Marquardt  parameter  (X)  DDAM,  EIGMIN 
computes  the  eigenvalues  EV(4)  and  eigenvectors  EGV(4,  4),  and  the  outer 
product  EH(4,  4,  4),  for  J.  EIGMIN  then  computes  X^  such  that  the 
Marquardt  matrix  is  positive  definite,  and  the  correspond¬ 
ing  parameter  correction  factors  C0RR(4)  =  +  X^C^j  g^. 


8.  Subroutine  TRIDMX  (N,  NM,  A.  D,  B) 

Reference:  UNIVAC  Large-Scale  System  MATH-PACK,  Programmer’s 
Reference,  UP-7542,  Rev.  1,  Sec.  9,  p.  1. 

TRIDMX  transforms  a  real  symmetric  matrix,  B(4,  4),  to  tri¬ 
diagonal  form  using  Householder's  method,  where  D(4)  are  the  resulting 
diagonal  elements  and  B(4)  are  the  off-diagonal  elements.  Input 
integers  N  and  NM  are  equal  to  the  order,  4,  of  B. 


I 

r 
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Reference:  UNIVAC  Large-Scale  Systems  MATH-PACK,  Programmer's 


Reference,  UP-7542,  Rev.  1,  Sec.  9,  p.  8. 

EIGVAL  evaluates  the  eigenvalues  of  a  symmetric  tridiagonal 
matrix,  using  Sturm  sequences.  A(4)  are  the  diagonal  elements  and  B(4) 
are  the  off-diagonal  elements  of  the  matrix.  The  eigenvalues  E(4)  are 
stored  in  descending  order  of  absolute  value. 
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10.  Subroutine  EIGVEC  (LP.  NM.  R.  A.  B.  E,  V.  P.  Q) 


Reference:  UNIVAC  Large-Scale  Systems  MATH-PACK,  Programmer's 
Reference,  UP-7542,  Rev.  1,  Sec.  9,  p.  15. 

EIGVEC  evaluates  the  eigenvectors  of  a  real  symmetric  tridiagonal 
matrix  using  Wilkinson's  method.  A(4)  are  the  diagonal  elements  and 
B(4)  are  the  off-diagonal  elements  of  the  matrix.  E(4)  are  the  eigen¬ 
values,  and  V(4,4)  are  the  eigenvectors. 


11.  Subroutine  MARQ  (EEH,  EEV,  DDLAM,  CCORR,  GGGJ,  HHH) 

Given  the  outer  products  EEH(4,  4,  4)  of  the  eigenvectors  of  the 
Hessian,  the  eigenvalues  EEV(4),  the  gradient  GGGJ(4),  and  Marquardt 
parameter  (X^)  DDLAM,  MARQ  computes  the  step  CC0RR(4),  +  X^C^  g^, 

in  parameter  space. 
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APPENDIX  B 


RELY  I  GLOSSARY  AND  INDEX 

(Library  subroutines  RANDEX,  RANDU,  RANDN,  TRIDMX,  EIGVAL,  EIGVEC 
are  not  included  here.  See  UNIVAC  MATH-PACK  references  given  in 
Appendix  A  for  detailed  information.) 


MAIN  (including  internal  subroutines  EIGMIN  and  MARQ) 


Variable 

ALPH 

B(4,  4) 


Description  Line  Number 

18,  22,  25,  30 

Value  for  a  in  simulation  mode 

162,  167,  170,  172 
Normalized  Hessian  matrix  in  EIGMIN 


BETA 


18,  22,  25,  30 

Value  for  B  in  simulation  mode 


CNC 

CND 

CORR(4) 

(MARQ: 

COST 


24,  33,  34 

Cumulative  values  for  n  in  simulation  mode 

c 

23,  32,  34 

Cumulative  values  for  n,  in  simulation  mode 

a 

3,  48,  57,  58,  61-63,  72,  78,  86-88, 


95, 

114,  115-117, 

127, 

137, 

160, 

163, 

201 

Vector  of 

corrections  to 

parameters 

CCORR) 

205 

,  207,  209,  229 

25, 

30,  43,  46,  49 

,  50, 

56, 

65, 

67, 

69, 

70,  73,  76,  91 

,  93, 

94, 

100, 

120, 

122 

,  124-127,  129, 

130, 

135, 

,  145 

,  147 

Most  recently  computed  value  for  the  mean-squared  error 


COSTI  46,  130,  145 

Cost  for  currently  accepted  parameter  values  used  in  mode  B 
to  determine  whether  the  new  estimate  for  the  parameters 
reduces  the  cost 


DDD  194,  195 

Denominator  used  to  normalize  the  matrix  (H  +  XI)  in 
EIGMIN 
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Variable 

DDDD 

DIA(4) 

DLAM 

(in  EIGMIN 
(in  MARQ: 

DND(300) 
EGV(4,  4) 
EH (4,  4,  4; 

(in  MARQ: 
ENC(300) 

ERR 

EV  (4) 

(in  MARQ: 
GJ  (4) 

(in  EIGMIN 
(in  MARQ: 

H(4,  4) 

(in  EIGMIN 
(in  MARQ: 


Description  Line  Number 

223,  224 

Denominator  used  to  normalize  the  matrix  (H  +  XI)  in 
MARQ 


162,  170-172 

Diagonal  entries  of  the  tridiagonalized  Hessian  matrix 
used  in  EIGMIN  to  compute  eigenvalues 

45,  48,  72,  74,  78,  95,  103,  111, 
114,  128,  137 

160,  188,  205,  217) 

182) 


4,  25,  32,  43, 

67, 

76,  91,  120,  135 

Incremental  values  in  n, 

a 

162,  172,  176 

Eigenvectors  for  the  Hessian  matrix 

3,  48,  72,  78, 

95, 

114,  128, 

137, 

160,  162,  176, 

188 

Outer  products 

of  the  eigenvectors  for  the  Hessian 

matrix 

used  to  compute 

(H  +  XI)~1 

EEH 

205,  207,  217) 

4,  25,  33,  43, 

67, 

76,  91,  120,  135 

Incremental  values  for  n 

c 

5,  7,  56,  93, 

126 

Value  for  termination  criterion 

3,  48,  72,  78, 

95, 

114,  128, 

137, 

160,  162,  172, 

181, 

182,  185, 

188 

Eigenvalues  for 

the  Hessian  matrix 

EEV 

205,  207,  217) 

3,  25,  43,  48, 

67, 

72,  76,  91 

,  95, 

114,  120,  127, 

135, 

137 

Gradient  vector 

for  the  cost  functional 

GGJ 

160,  163,  201) 

GGGJ 

205,  207,  229) 

3,  25,  43,  48, 

67, 

72,  76,  78 

,  91, 

95,  114,  120, 

127, 

135,  137 

Hessian  matrix 

for  the  cost  functional 

HH 

160,  162,  167, 

194) 

HHH 

205,  207,  223) 

Value  for  X 
DDAM 
DDLAM 
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( 


l 


Variable 

I 

ICLAM 

I  FLAG 

INDEX 

ISIM 

J 

JDEX 

KK 

KL 

LMBEX 

NA 

NJ 

NN 


Description 


Line  Number 


Index  for  various  loops 

104,  113,  144 

Index  which  counts  the  number  of  times  that  X  is  increased 
in  mode  B 


47,  70,  71,  106 

Index  that  counts  the  number  of  times  that  an  iteration  of 
mode  A  reduces  the  cost  by  less  than  1  percent 

51,  80,  81 

Index  that  counts  the  number  of  times  that  the  step  size 
has  been  reduced 


15,  17 

Indicates  whether  the  run  is  a  simulation  (ISIM  =  1)  or  an 
estimation  with  real  data 


Index  for  various  loops 

105,  143,  144 

JDEX  *  1  indicates  the  first  time  that  changing  X  has  been 
successful  in  a  given  iteration  of  mode  B 

Index  used  for  DO  loop  in  EIGMIN  for  computing  outer 
products  of  eigenvectors 

184,  185,  188,  214,  217 

Index  used  for  DO  loop  in  EIGMIN  and  MARQ  for  computing 
(H  +  XI)-1 

112,  122,  123 

LMBEX  =  1  indicates  that  X  was  changed  in  mode  B 

2,  18,  22,  25,  30 

Value  for  N  in  simulation  mode 
a 

10,  25,  43,  67,  76,  91,  120,  135 
Number  of  test  intervals 

5,  7,  10,  11,  22,  25,  31,  43,  67,  76, 
91,  120,  135 
Number  of  tape  versions 


Variable 


Line  Number 


Description 


OFDI (4) 

162,  170-172 

Off-diagonal  entries  in  tridiagonalized  Hessian  matrix 

PCOST 

50,  69,  70,  73,  9A,  12A,  127,  129 
Current  minimum  value  for  the  cost  functional 

PHI 

18,  22,  25,  30 

Value  of  <S  in  simulation  mode 
a 

R(4,  A) 

163,  166,  188,  195,  201 

Matrix  (H  +  XI)  computed  in  EIGMIN 

RR(4,  A) 

207,  211,  217,  22A,  229 

Matrix  (H  +  XI)  computed  in  MARQ 

REFF 

26,  30,  AA,  A9,  68,  77,  92,  100, 
121,  122,  125,  136,  1A7 

Estimated  number  of  errors  remaining  at  the  end  of  test 
period 

RMTTF 

25,  30,  A3,  49,  67,  76,  91,  100, 
120,  122,  125,  135,  147 

Estimated  mean  time  to  failure 

RRR 

18,  19,  22 

Randomization  value  in  simulation  mode 

S  (300) 

3,  12,  22,  25,  38,  43,  67,  76,  91, 
120,  135 

Error  data 

T (300) 

3,  13,  25,  43,  67,  76,  91,  120,  135 
Tape  version  replacement  times 

TAU 

57-63,  79,  86-88,  102,  115-117,  144 

Step  size 

TD 

5,  7,  13,  22,  25,  43,  67,  76,  91, 
120,  135 

Length  of  each  test  interval 

TEMPl(A) 

TEKP2(4) 

163,  171,  172 

163,  171,  172 

Vectors  which  are  used  temporarily  in  the  computation  of 
the  eigenvalues  and  eigenvectors  of  Hessian  matrix 
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Variable 


Description 


Line  Number 


ZA 

ZB 

ZN 

ZNC 

ZP 

ZZA 

ZZB 

ZZN 

ZZP 

A  (dbl) 

B  (dbl) 
D 


40,  43,  49,  52,  67,  82,  91,  96,  100, 
107,  120,  122,  125,  131,  135,  139, 
147 

Value  for  a  in  simulation  and  estimation  mode 


40, 

43, 

49,  53 

,  57, 

,  61, 

64, 

66-67, 

83, 

86, 

89-91, 

97, 

100, 

108, 

115, 

118- 

-120, 

,  122, 

125, 

132, 

135, 

140, 

147 

Value  for  8  in  simulation  and  estimation  mode 

40,  43,  49,  55,  63,  67,  85,  88,  91, 

99,  100,  110,  117,  120,  122,  125,  134, 
135,  142,  147 

Value  for  N  in  simulation  and  estimation  mode 
a 

26,  30,  44,  49,  68,  77,  92,  100,  121, 
122,  125,  136,  147 

Estimated  value  for  nc  at  the  end  of  the  test  period 

40,  43,  49,  54,  58,  62,  67,  84,  87, 

91,  98,  100,  109,  116,  120,  122, 

125,  133,  135,  141,  147 

Value  for  <t>  in  simulation  and  estimation  mode 
a 


52, 

76, 

82, 

96, 

107, 

137, 

139 

Currently  accepted 

value 

for 

a 

53, 

76, 

83, 

97, 

108, 

115, 

132, 

140 

Currently  accepted 

value 

for 

8 

55, 

76, 

85, 

99, 

110, 

117, 

134, 

142 

Currently  accepted 

value 

for 

N 

a 

54, 

76, 

84, 

98, 

109, 

116, 

133, 

141 

Currently  accepted 

value 

for 

*a 

Subroutine  RNDTA 
1,  7,  19 

Correction  rate  parameter  a 

1,  8,  19 

Correction  rate  parameter  8 

4,  21,  18 

Estimated  number  of  remaining  errors  nfl 
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Variable  Description  Line  Number 

E  (dbl)  17,  18 

D 

JJJ  1,  15 

Number  L  of  intervals  (software  versions) 

15,  21,  27 

Index  corresponding  to  k.th  interval 

1,  3,  10,  16,  17,  19 

Initial  number  (y-augmented)  of  software  errors 

1,  9,  16,  19 

Detection  rate  parameter  cf> 

a 

4,  7 
A 

RB  4,  8 

B 

RNA  4,  10 

NA 

RP  4,  9,  21 

P 

RQ  5,  6,  21 

RR 

RR  1,6 

Storage  place  for  MAIN  input  initial  random  number  for 
RANDN.  Contains  first  exponential  random  number  from 
RANDEX  upon  return. 

RT  4,  11,  21 

T 

4,  21,  22,  27 

Poisson  random  4^(10  returned  by  POISS 

1,  11,  16,  19 

Time  interval  T 

V(300)  (dbl)  1,  12,  27 

RZ,  Poisson  random  sequence  A^Oc)  returned  by  RNDTA 


RZ 

T  (dbl) 


K 

NA  (dbl) 

P  (dbl) 

RA 
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Variable  Description 


Line  Number 


X(2) 

(dbl) 

Temporary  memory 

for 

12-14,  16,  17,  19,  20,  22 
current  cumulative  random  N 

25 

Y(2) 

(dbl) 

Temporary  memory 

for 

12,  16,  19,  20 
current  cumulative  estimate 

-► 

M 

Subroutine  POISS 


C(100) 


DD 

K 

KKK 

PP 

Q 


RRRR 


TT 

U 

ZZ 


2,  4,  6,  8,  11 

Set  of  uniformly  distributed  random  numbers  generaged  by 
RANDEX  to  be  used  by  POISS  as  the  sequence  of  times 
between  successive  detected  errors 

1,  7 

Number  of  remaining  errors  n^ 

10,  11,  13 

Counter  of  successive  detected  errors 

1,  3 

Number  L  of  test  intervals 


Parameter  value 


1,  7 


9,  11,  12 

Cumulative  time  I^C^  during  the  test  interval,  accumulated 
until  it  exceeds  T 


1,  * 

Initial  random  number  for  starting  RANDN.  Its  value  for 
k  =  1  is  input  by  MAIN,  Subsequent  values  are  set  by  POISS, 
RRRR  »  225  c(100). 

1,  12 

Test  interval  T 

7,  8 

Mean  frequency  <t>ana  t*ie  error  detection  r^ 

1,  13,  16 

Poisson  random  number  of  detections  A<j(k)  generated  by 
POISS  for  the  k^  interval.  It  is  Aj  such  that: 


l  C  <  T  ,  A  <  100 
i-1  1 


* 

4 

i 


i 
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APPENDIX  C 


PROCEDURE  FOR  OPERATION  OF  RELY  I 


The  program  listing  (FORTRAN  V,  Appendix  D)  accompanying  this 
report  consists  of  MAIN  with  internal  subroutines  EIGMItJ  and  MARQ,  and 
external  double-precision  subroutines  RNDTA,  POISS,  KOST,  TRIDMX, 

EIGVAL,  and  EIGVEC.  The  subroutines  which  call  single-precision  library 
functions  have  the  necessary  coding  for  converting  between  double-  and 
single-precision  variables.  Certain  double-precision  library  functions 
are  used,  viz.,  DABS,  DEXP,  and  DSQRT. 

The  input  deck  depends  on  whether  data  are  to  be  simulated  or 
are  to  be  read  from  input  cards. 


INPUT  DECK  (to  simulate  data  and  estimate) 

Card  1:  TD,  NN,  ERR,  [F5.2,  2X,  14,  2X,  F8.6] 

Card  2:  ISIM  [12]  (must  be  unity) 

Card  3:  ALPH,  BETA,  PHI,  NA  [4(G14.6,  2X) ] 

Card  4:  RRR  [16X,  G14.1]  (input  initial  random  number) 
Card  5:  ZA,  ZB,  ZP,  ZN  [4(G14.6,  2X)] 


INPUT  DECK  (to  estimate  using  punched  card  data) 

Card  1:  TD,  NN,  ERR,  [F5.2,  2X,  14,  2X,  F8.6] 

Card  2:  ISIM  [12]  (must  be  zero) 

Card  3:  ZA,  ZB,  ZP,  ZN  [4(G14.6,  2X)] 

Card  4-(k+3):  S(i),  i  -  1,  2,  ...  K  [G10.1] 

The  integer  part  of  any  real  number  RRR,  0  <  RRR  <  2^f 
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determines  a  unique  repeatable  random  sequence  A, (k)  from  the  simulator. 

d 

-► 

Initial  guesses  6q  for  the  model  parameters  are  recommended  as 
follows: 

0  <  ZA  <  1.  (typically  0.8) 

0  <  ZB  <  1.  (typically  0.0) 

0  <  ZP  <  0.2 

0  <  ZN 

To  produce  different  simulated  random  data,  the  operator  must 

35 

change  the  input  initial  random  number  RRR,  0  <  RRR  <  2 

Though  J  <  ERR  is  the  internal  stopping  criterion,  experience  in 
random  cases  proved  maximum  pages  of  printed  output  (say,  10)  to  be  as 
practical  a  stopping  criterion  as  any. 

The  program  prints  out  the  current  accepted  parameter  estimates 

*► 

0  “a,  6,  4>,  N  ,  together  with  running  estimates  of  the  reliability 

parameters  NR  and  MTTF,  and  selected  auxiliary  quantities,  at  each 
step  i.  However,  the  label  "CHANGING  LAMBDA"  indicates  only  tentative 
parameter  values  produced  in  mode  B  (see  Sec.  3.2).  Therefore,  these 
tentative  values  must  not  be  taken  as  values  which  minimize  the  cost 
functional  J.  The  final  estimates  of  the  reliability  parameters  are 
those  associated  with  the  last  accepted  iteration  i. 
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APPENDIX  D 


HtL  Y*ht.l.  Y  l  i  ) 
1 
2 
2 
4 
b 
o 
7 
o 
9 
10 
11 
12 
12 
14 
13 
lo 
17 
lo 

19 
tU 
21 
<12 
cl 

it 

2b 

20 
27 
20 
29 

20 

21 
22 
22 

24 

25 
2t> 

27 

2b 

29 

40 

41 

42 
42 

44 

45 

4b 

47 

4b 

49 

50 

51 

52 
52 

54 

55 
5b 


RELY  I  PROGRAM  LISTING  AND  SAMPLE  OUTPUT 


•  HA  14 

l.-PLICI  I  REAl*8(a-H.0-Z> 

KEAL*8  i*A 

U1MENSIU.  T (  200  )  ,  S  <  200  >  »h  {4 . 4  >  »t»J(  4  )  ,C0RR  ( n  )  ,EH  <  4 ,4 ,4  >  ,EV(4) 
blMEMSIuN  QNL  <  3Uo  1  >  tNC  ( 2vu  1 
KEAu  ( 5»  ob  7)  TL.Nh.LRR 
bu7  FoRMAT(Fb.2#2X. l4.2X.E8.o) 
wk1Tl16»836)Yo  i  HK  .ERR 

006  F0RRAT ( jKi • Interval  LENGT n  =* .2A.F5.2.2X. ‘NUPPER  oP  INTERVALS 
1cX*14.2X.  •TERKINmUON  CRIIEkION  =«.FR.6) 

No=NN 

UO  14  J— 1 . NN 
S(J)=0 

l  (J)  =  (J“i)*TLi 
14  CONTINUE. 

KtAblb.ibO) IS1M 
150  KuRklAK  12 1 

if  (ISIK.hE.UGO  10  10 
«c.AO  ( b»cfa8 )  AcPH.olTA  »FHI  •  uA 
nLAU ( 5.oe9 ) KKR 
*hITt (6»b89)RKR 
ae9  KuKKATUbX. 014.4) 

CALL  RNoTAIAlPH.cETA.PHI . „A » TD . S . P°R . MN ) 

C..0=0 

C  <C=0 

CALL  K  Obi  ( ALPH » fatT  A  »PHI »  i»«  » S.NJ*  T  .NN.TD.H.  GJ»  PMO  »  COST » FllC  •  °MTTF . 
l2i«C  »kERk ) 
t>R  I TE  ( 6  > i69 ) 

lo9  F  ORHAT  ( 27X  i  •  xLPhA  • » 5X  * •  tii-TA  »»4X,»  PHI  *»5X.»  NA  '.5X, 

1*  COST ,»2X»,CST.NC,»4X»'  WTTF  * , 4X . »I\IA-NC • ) 

*KlTE(6.170)  ALPh.uETA.PHI  .NA  .COST.  7NC  .RMTF.RERR 

L'O  89  1  =  1  .NN 

C,,u=CNC+LNO(I) 

U.C=CNC+LNC<  I) 

890  )  CNO. LNC 

aVO  F ORF  AT  (nXi  ' NG=  •  » cX ,F  10. 2 .  22 • ' NC=  * . 2X . FlO . 3 ) 
o9  CoNTIf.Ut 
uu  TO  ll 

10  RLAU<5. 151) (S( I) .1=1. NN> 

Ibl  FORMAT ( ul 0 • 1 ) 

11  RLAU ( 5.086 )ZA» Zt)»2P*ZN 

886  F'GRF,AT<bl4.6»2X«bl4.6»2X»ljl4.6.2X.G14.6> 
akITE (6* 169) 

CALL  KObTIZA.ZB.ZP.ZN.S.Nw.T.NN.TC.H.GJ.DHP.COST.ENC.RRTTF. 
lZNC.RERh) 
uLAK=l 
COSTlsCcST 
IFLAOSO 

CALL  EIuMIN IH.COKh » GJ.DLaM.EH.EV ) 

72  taRlTt(6’171)ZA.2b.ZP.ZN.CUST.ZNC.R»TTF.RERP 
75  PC05T=COST 
lHOtX=0 
ZZA-IA 
225=26 
ZZP-ZP 
Z2N=ZN 

IFICOST.LT, EKRIOO  to  73 


37 

30 

39 

00 

ol 

oi 

63 

04 

ob 

00 

o7 

0O 

04 

7u 

71 

7c 

73 

74 

75 
7o 
77 
7o 
79 
oO 
01 
0C 
03 
04 

eb 

00 

07 

0b 

09 

9U 

41 

9C 

43 

94 

4b 

90 

97 

90 

99 

lub 

lbl 

1"* 

Hi 

104 

105 
100 
107 
loo 

109 

110 
111 
1U 
Hi 


TaU=  (5-4.0 )  / (2*DAfcb (CORR  (2)  >  ) 

1F(  I AU.ol .  ( .C-ZP)/(2*0AbS(CoRR(3)  >  )  >  TAU=  (  .2-ZP)  /  <2*DABS  (COPR  ( 3) )  ) 

if-  <TAU.0T.1JTAU=1 

IF(TaU. 01. .0001)00  TO  74 

Zb=2B-CuKR(Z)*TAU 

Ap=2P-CuhR(3)*TAU 

Zi.=ZN-C0Krt(4)*TAU 

IF (Zb. LI .0)Zb=0. 

1F(C0ST.GT.1Q)0LAF,=  1. 

IFUb.GT  .10.)ZB=1G. 

OAU.  ACsT  (ZA»ZB#Ap.ZN#S.Nu.T.tiN»Tr'.H»GJ.DNP.COST.EFIC.RMTTF. 
iZf.C.REKK) 

1F(C0ST.GE.PC0ST)G0  TO  70 
IF (COST/PCOST.GT. .99) IFLa3=1FLAG+1 
IF ( lFLAo.GT .4)  Gu  TO  74 
CALL  EI»MlN(Hi COAk  > G J •  DLAi*! » tH .  E  V ) 

KCOST^CobT 
DLAM=LLmK/10 
OO  TO  7c 

76  CALL  K03l  ( ZZA .  ZZb .  ZZP.  ZZU.S.Nj » T  . NN , TO .H. G J.OUD* COST »£NC .RMTTF* 
lZf.C.KtPn) 

C«LL  £  I<jF  IN  (H  •  COKK  #  GJ  i  DLAi"' »  0H  t  EV  ) 

71  TA’J=TAU/10 
Ii.GtX=J..uEX*l 
IF ( IM/Ea.GT .b)Gu  TO  74 
ZA=Z2A 
Zu=ZZb 
cP=ZZP 
ZUSiZU 

Zo=ZB-COkK<2)*Tao 
aP=ZP-CGKR(3)»TaU 
2i  tsZN-Cokrt  ( 4 )  *TaO 
IF'  (Zb. cl  .0)ZU=0. 

IF U0.GI. 10. 126=10. 

Call  KObT (ZA.ZB.ZFiZNiS.Nuil .DN.TD »H .GJrDMD .COST .ENC .H^TTF . 
lZ.MC.REHk) 

IF (COST .LT .ERR ) GO  TO  73 
IF  (COSI .GE.PCObT )G0  TO  7l 
CALL  £Ioia1N(H>C0HR.GJ>0LAm>£H>EV) 

ZZASZA 

cZB-Zb 

ZZP=ZP 

ii.US.ZU 

*hITt(6*172)ZA.Zd.ZP.ZN.C0ST.ZNCfRMTTF»RERP 
60  TO  7a 
74  TAU=.l 
ULAMS.Ol 
ICLamsO 
JUEXSO 
1FLA6S0 
ZA=ZZA 
Zb=ZZB 
ZP=ZZP 

iUslZn 

77  uCAM=10*OLAM 
LPBtXSl 
1CLAMSHLAM41 


11* 
113 
11*3 
117 
110 
11* 
ltU 
121 
122 
Lei 
12* 
l2b 
l2o 
1*7 
Leo 
12* 
IjO 
ljl 
1  o2 
133 
lb* 
133 
loo 
1  j7 
lou 
139 
1*0 
1*1 
l*c 
1*3 
1** 
1*5 
l*o 
1*7 
l*o 
1** 

150 

151 

152 

153 
13* 
135 

150 

137 

156 

139 

100 

101 
lot 
163 
lb* 
103 
loo 
107 
100 
lo9 
170 


C^UL  MAh*  (EH, EV  ,  ULAM » CORK  *  *3  J  *  H ) 

76  4o=2Zb-CuRR(2)*TAU 
2P=aZP-L0RR(3)*TA0 
ZI,=ZZN-C0RR(*)*TAU 
IFUb.GT. 10)26=10. 

IF (tb. LI .0)26=0. 

CALL  KOST (ZA»/B.2P«ZN»S»Nj»T »NN > TD ,H » 6 J.DND » COST ,£NC»RMTTF . 

12NC  >KERk) 

IF (LMDEa.EQ. 1) WHITE (6.190 )ZA. ZB. ZP.ZN.COST.ZNC.RMtTF.RERR 
LVBtX=0 

IF  (COST. GT.PCCS'i  1*0  TO  77 

*KlTL(6.191)ZA.2t,.2P,ZN.CoST.2NC.RMTTF.RERP 
IF(COST.LT.EHR)GC  TO  73 

IF  (COST /FC0ST.LT,,9.Af!0.Pi,03T.GT  .SO)  CALL  EIGMlH (H .CORK » GJ. 
1l,lAM.EH.lV) 

IF (COST/PCOST .LT . .9. ANO.PCOST ,6T ,50)GO  TO  72 
IF (CCSTi .GT.CGST ) GO  TO  91 

4A=4lA 

2u=220 

4M=ZZP 
4 U-ilU 

C„Li-  KObT  (ZA.2B.2P.ZfJ.S.r.„.T.NN.TD.H.GJ.DMn.C05T.LNC.KMTTF. 
12i.C.Htr<K) 

Call  CIoMN(H,COKR.GJ.OLAk.EH»Ev> 

0*  TO  7t 
*1  Z2A=ZA 
ZZb=ZB 
22P=ZP 
Z2N=2N 
JOEA=JCtA+l 

iF (wOEX.tQ.l)TAU=TAU»5**ICLAM 
C0Sll=Cu5T 
GO  TO  7b 

73  nRlTE(6.173)ZA.2b>ZP.ZN.CoST«ZNC .RMTTF . RERR 

170  FORMAT  (s*.  •SIt-iULATlON  VALvEs  ARE  * .  2X  ,  3  (Ffl .  5 . 2X )  , 

13(Fo.2.4X).Fltl.b.2X.Ffi  .2) 

171  format (ox . *Ntw  parameter  values*, ex, 3(Fb. s. 2X> » 

13(F6.2»4X) .F10.O.2X.F8.2) 

172  FORMAT (3X» ‘AFTER  REDUCING  STEP  SIZE* ,4X,3(F8.5»2X) , 

13(Fb.2.4X) .F10.0.2X.FB.2) 

173  FORMAT (*X, ‘FINAL  VALUES  AhE* »2X.3(F8.5,2X) , 

13(F6.2.4A) .F10.O.2X.F8.2) 

190  FoRMAT(3X, 'CHANGING  LAMBDA  • ,2X.3(F6.5.2X) . 

13(Fo.2.4X) .F10.6.2X.F8.2) 

191  F0RMAT(5X. ‘STEEPEST  DESCENT  VALUES  * »2X,3(Ffl.5.2X) » 
13(F6.2»&X)>F10.G*2X.F6.2) 

SUBROUT l NE  E IGM IN ( HH . CORR .SGJ.DOAM.EH.EV) 

IMPLICIT  REAL*8(A-H»0-Z) 

DIMENSION  HH ( * » * ) ,B(*<*> ,tV(*) ,EGV(4,*),EH(*»4,4) ,OIA(*> ,OFDI(*> » 
1 T  EMP 1(4). TEMP2 ( 4 ) ,R ( 4 .4 > .CORR (4 > ,GGJ (4 ) 

DO  700  1=1,4 
DO  701  w=l»4 
R ( I , J>  =U 

b(I«J)=MH(I,J) /DSLRT ( DABS (HH ( I , I ) *HH ( J. J) ) ) 

701  CONTINUE 
700  LONTINUt 

CALL  TRlUMX<4»4,b«DIA.0FDI) 


1?1 

CALL  EIoVALl4»E\/.LIA.0FCI.TtMPl.TEI'P2) 

17* 

CAUL  EI«VEC(4,4.b.0IA.0PDi.tV.EGV,TEMPl.TEMP2> 

174 

DO  708  1=1.4 

174 

DO  709  d=1.4 

17b 

uO  71C  kN=1<4 

1/b 

EHU.J,K*.)=EG\,  (u.I)»EGV(KK,i) 

177 

710 

COliTIf.Ut 

176 

7o9 

CONTINUt 

17V 

/ofl 

CONTINUt 

loO 

DC  713  1=1.4 

lei 

IF(EV<I).GT.O)GO  10  713 

16* 

IF (DDAM.lT . “LV ( 1 ) ) D0AP=DAuS (EV ( i ) ) l . 1 

loo 

713 

Continue 

lo4 

DO  70*  Nt=l,4 

lob 

1F<dm&SIUV(KL) ).LT..0001)tV(KL)=l 

loo 

DU  703  1=1.4 

107 

bC  704  b=l,4 

lob 

MI.  J)=K(I.UJ*Eh(KL»  I.  Ji/lEoCKUM-DDAMi 

le9 

704 

CuNT INUt 

190 

7o3 

CoNTINUt 

191 

7l)2 

CONT  IMUc. 

IV* 

ul  705  i=1.4 

19J 

DC  706  o=1.4 

194 

ojO=USGkT (OAbSCHrU  I . I )*HH< J. J) ) > 

19b 

MI,  J)=MI.J)/DuD 

19o 

706 

CoI.rK.UL 

197 

705 

CONTIliUt 

19b 

UO  711  1=1.4 

199 

CoRK ( I ) =0 

*00 

Oo  712  v=1.4 

201 

CORK (I) =CORR ( I ) +K ( I . J ) *6Go ( o) 

*02 

712 

continue. 

200 

711 

CONTIliUt 

*04 

kETcRN 

*0b 

budKOUT ll.E  MAKQ (tEH . EEV.DDUmM.CCORP . GG6J.MHH ) 

200 

IMPLICIT  REAL*8(A-H.0-Z) 

*07 

DIMENSION  EEH  (4*4.4)  . EEV  ( 4  i  . CCORR  ( 4 ) « G(»GJ  ( 4 )  »RM  (4.4)  .HHH  (4,4) 

*0b 

UO  714  1=1.4 

209 

CCOKK ( I ) =0 

210 

DO  725  *1=1.4 

211 

KH(1.U)=0 

212 

7*5 

CONTIliUt 

214 

714 

CONTINUE 

214 

DO  71b  r>L=l»4 

21b 

UO  717  1=1.4 

21o 

UO  716  d=1 *4 

*17 

RK l I » J) =KH l I . J) ♦tEHIKL. I  ,*,)  /  (EEV (KL) +COLAM) 

210 

718 

CONTIMUt 

219 

717 

CONTINUt 

220 

716 

CONTINUt 

2*1 

UO  719  1=1.4 

2** 

UO  720  J=1.4 

2*4 

UUDU=CS4nT (UABS (hHH ( I . X ) *nHH ( J» J) » ) 

2*4 

KR( I , J)=HR ( I ,U)/UUDO 

2ib 

720 

CONTIliUt 

2*o 

719 

CONTINUt 

2*7 

UO  721  1=1,4 

! 


2tb  Lu  722  w—  l » 4 

CCOf<R<I)=CCOkR<X  )»Rh<I»J).S6GJ(J> 
2o0  722  CCKUMJc 

Hi  721  tONTIHUt 

2*2  k£Tu«r. 

2-i *  fc.NO 

<»•> 


i»khT*S  r. tk-J  ( 1  i  ,/\OST 
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fit 4.Y  ««£ O  l  ± )  ,KoST 


2 

4 

4 

3 

o 

7 

o 

9 

1U 

11 

12 

14 

14 

x’j 

lo 

n 

lu 

j.* 

cu 

cl 

kd 

c4 

24 

CO 

2o 

d7 

CO 

c9 

44 

4i 

41 

44 
04 

45 
4o 
47 
oo 
49 

HO 

Hi 

42 
Hj 

44 

46 
Hb 

47 
46 
49 
60 
61 
62 

64 
34 

65 

6b 


ioBKOt'Ti  I»£  KOS  7  C  A  A .  83  •  PP  .  i.N*  .  SS  .  fl  J  J  *  TT  ,  NMN  .  DT  *  H .  G  J  •  UNO  • 
lCoSr#CN4.R«TTFi2NC»REPR) 
implicit  keal»8u-h.o-2) 

KcAL«B  I’.l'iA 

UIMtliSIoN  TT<300  .55(300)  ,Go(4)  *H(4,4) 
t>lKfc.jSIwi\  0Um0<4)  i  ONQ  ( 30Q ) « uHC  (400 1  .z*'TTF<  jooi 
uIMcnSIoN  01(4.4) 

CoST=0 

ZLLUSl 

cLc2<51 

24.02=0 

ccL2l=0 

2!.L=0 

cl.CsO 

oO  14  .1=1.4 
oLf.u(I)  =  0 

oJ(l)=0 

uC.  12  J=1.4 

H l I «  J) Su 

LI ( i . J)=U 

12  bGUT  If 'bt. 

13  LUN  r InUt 
Ml=-AA*f-r'*QT 

Ac=<  1-Am)  «CT*b04,.A-PP*0T*AA**2 
2L=l-0T*laB*AA*HK) 

«7=PP*Uf 

•*«J=  (FP*0T)  **c 

A9=AA *Mb 

A1u=A6*mA«*2 

AU=AA*a7 

Al^sFP*  ( AA*0T )  **.. 

bl=t4e*oT4AA-AA*Ol«(ar!4AA»PP) 

UJNu(3)=l«NA*jT 

uuNu(4)=pP*0T 

04  10  K  =  l.f\itoN 

o.\Q  »* )  sr-P*Nf«A»0r*  <  ZUL1  !♦  A**ZLL12 ) 

UlCIM  =FP*NNA*DT* (  ZU.21+A**«:LL22 ) 

2r<u=2ND*bN£JtK) 

/:«C=ZNC-*bNC(K) 

Z.iTTF  (K )  si/  (PP*  (NNA-ZNC  I I 
1F(K.E«.NJ<N)R*TTF=Z*TTF(K) 

IF  (K  .EG.r.HN)  RERR=f;NA-ZNC 
1 F  ( K .  L 0  .  i«MN )  ZZNASUN A -ZN C 
COSTsCObT* (SS (K )“1NU (K ) ) »*2 
JO  41  1=1.4 

6 J  ( 1 )  «“2*  ( SS  ( K  1  -ul<0  ( K  )  1  •UuNO  ( I  >  FG  J  ( J ) 

JO  42  L=1 » 4 

M(  I*O=2*(DDND(l)*0JND(L)-(SS(K)-0N0(K)  )*01  <  I.U  UHd.L) 
42  CONTI NUL 
41  CONTINUE 

ZLKll=Z4Lll*ZLL12*Bb*0T 

2l.K12sZi.Ll  1*  ( -PP*lT  >  *ZLL1^«2L 

2LK21sZLL21*ZLL22»BB*0T 

2LK*2=2u.21*(>PP*UT)*ZLL22*2L 

oLNJ(1)=PP*NNA*DT*(K*A1*ZLH2*ZLK12) 

ooNJ(2)=PP*NNA*OT*(K*A2*ZlH2) 
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b? 

LUiOl  3)=nNA*oT*(ZuKU4AA»2LM2)+PP*NNA*0T*K*  (-A  A)  *0T* 

bo 

1(<.L011+«*.*ZLL12) 

b9 

uoN0!4)=PP«[jT*(ZLMl-*AA*ZoK12) 

o(J 

lF(K.r,£.l)G0  TO  3u 

ol 

CJ 1  ( 1 .  3 )  =-2*PP»NUA*0T 

02 

Ul ( i » 4 ) =-UT*PP»«2 

oj 

01(3*3) =-2*AA*Nh«*GT 

o4 

01 (3»4)=UT-Z»AA«PP*UT 

ob 

ul(3*l) =01(1*3) 

oo 

01(4*1) -01(1*4) 

o7 

U 1 ( 4  *  3 ) — 01(3*4) 

oo 

00  TO  31 

o9 

30 

01 (1*1)=PP«NNA*QT«K*( (K-l)*iyi2*A9-2*ZLL12*A7) 

7u 

01  (1*2)a*  P*(*mA*oT*K*(  (K-U*ZM12«A12-DT*2LL12) 

71 

ol  ( l*3)=i.!*A*0T»(K*ZLL12*(-All)i2LK12)+PP*MNA*0T«K*(  (K— 1 )  • 

72 

1  ( >>ii  1  *  A± i *DT*ZM12* AID  )  _2*tl.ol2*AA*r)T-ZLLl  1*DT  > 

70 

01(1*4) =PP*0T* ( a»2LL12* ( “All ) +ZoKl ? ) 

74 

04  (2*2)=PP*Nf4A*UT*K*(K-l)*ZM12*ijl»(“PT> 

7  b 

ol <2*3)=iw!A*oT*A»2LL12*bl4PF*NNA*DT*K*(K-l)*(ZMll»<-|jT>* 

7o 

1u1-*>12«aA«LT»B1) 

77 

u1(2.4)=FP«UT»K«4LL12*B1 

la 

«j1(3.3)=4*W<A*DT*A*(-AA)*uT*(  ZLlI  1  +  AA*ZLL12)  ♦  PP*N|*A* 

*9 

1  ( AA»«2)  *  luT»*3)  »l'.»(K-l  )*(4yil  +  AA*Z*'12) 

Oli 

Ll( J.4)ioT»(ZLKll+2LK12«4M)*PP*oT*K»(-AA)*r!T*(ZoLlJ4AA*ZLL12) 

ol 

kji(4»4)— 0 

o2 

01(2*1) =U 1(1*2) 

Ob 

ol(3*l)-o 1(1*3) 

o4 

01(4*1) — U 1(1*4) 

Oil 

wl ( 3  *2 ) Pul (2*3) 

OO 

01(4*21=01(2*4) 

O  7 

ul(4«3) =01(3*4) 

od 

01 

Zi-ll=ZLoll 

09 

zr..l*.=ZLol2 

90 

a.-i21=ZLo21 

91 

4-22=ZLL42 

92 

4LL11=Z4-K11 

93 

ZL012=ZOhl2 

94 

aLG21=Zo(\21 

9b 

aL022=ZoK22 

9o 

10 

CoNT INUt 

97 

ketusn 

90 

C...0 

<••> 

fcPKT.b  kE*.Y(1).SUd2 
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RtuY*i<ttVu)  .hub2 
1 
2 
i 
4 
3 
o 
7 
a 
9 
10 
11 
12 

13 

14 
lb 
lu 
17 
lb 
1* 

<u 

*1 

2b  08 

to 

c.7  UiO 

2d  o7 

29 

30 


300KOUT ll.£  Ki<LT A(A*titPtNAiTtVrRR* JJJ ) 
l.\PLICn  R£*L*4U-H.0-Z> 

K£AL«b  NA 

KLAL  KA * HS t KP $ RNX • RT iDt K2 

hEAL  ro 

RO=SNGL(h«) 

hA=SNGLl».) 

Kb=SN6L<b> 

KPbSUGLlH) 
hi  ,A=SNGl.  ( NA ) 
nT=SUGLlT) 

*RlTt(6rb7) 

wlMENSIbN  X(2) • V IbOO) *  t  (2) 

A(1)=0. 

X  1 2 ) -0 , 

UO  100  K=l<JwJ 
1  ill=X(i) ♦P»T»(itf-X(21) 
t=NA-X(t) 

3=SNGL(fcJ 

Y <2)=B*1*X(1)+(1-T«<L*A*P> )*X(2)+A*P*T*UA 
a(2)=YU) 

Call  POaSS (L) > R2 * hh > AT •  RG* A ) 

A  ( 1)=X ( i ) tOoLL (ht) 
i,f.  ITE  1 6  *  ba )  X  1 1 )  >  a  (2 ) 

roRMAT ( ±0X» *hU’ #2A.G14,6.*X*  *RC* .2X.G14.6) 
w(K)=06Lt(RZ) 

COUT l  MIL 

FORMAT^.  *RaN0OM  VALUES  FOk  NO  AND  NC  GENFRaTED  FOR  SIMULATION*) 

hLTuHN 

tt.U 
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CPkT.S  htuY (1) .SUbd 
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RLi.  Y  *KcL  Y  l  a  1 

l  • 5ub3 

1 

SUbKOUT li.E  RuISS  <  LO » Z2  *PP  >  TT » RRRR » KKK  ) 

£ 

OIMLnSIoi.  C(luO) 

O 

IF  (KKk.GT.DSO  TO  555 

4 

C  U)=ftRKR 

3 

60  To  6oo 

0 

355 

C ( 1 )=C ( aUO) *£*«£b 

7 

oofa 

O=00»FP 

o 

CALL  R AxoEX ( C  > 1 00  >  U ) 

y 

030  ■ 

10 

LO  100  n=1>1o0 

il 

0=G+C<K) 

12 

IF  (O.LT. TTIGO  TO  100 

13 

£Z=FLUAi IK-1) 

14 

60  TO  2uu 

13 

luO 

OONTIliUc. 

16 

22=100 

17 

<.00 

rtLTuRf* 

lo 


t,.b 


TRIOOCIO 


«T  ( 1 ) 

1 

2 

3  c 

4  C 

5  C 
b 

7  C 

a  C 

9  C 

1U 

n  lu 

12  c 

13  C 

14  C 
lb 

16  C 

17  C 

16  lb 

20  C 

<1  c 

22  C 

23  C 

24 

2b 

26  20 

27  C 

26  C 

29  C 

30  C 

31 

32  C 

33  C 

34 

3b  C 

3e  C 

37  C 

36  C 

39  C 

40  24 

41 

42 

43  3w 

44 

4b  37 

46  C 

47  4 

46  C 

49  C 

bo  C 

61  V 

62  C 
bo 

b4 

bb 

bo  o4l 


Souo 

SUBROUTINE  ThIOMA  (N.NH.fc.D.O) 
IMPLICIT  K£AL*ai*-H»0-2) 


TRluIAGONALIiATIGN  OF  REAL  SYMMETRIC  MATRIX. 


UiMENSIuN  AtWM.NMk  ,0(NMI  .uINMl  TR1D0020 


SAVE  ORIGINS.  UlAGON/LS  IN  ARRAY  D. 


UO  10  Isi.N  TRIP0040 

UU)=A(l.I)  TR 100050 


FOR  N-2  RETURN  VlTnOUT  COMPUTING. 


IF(n-2)ou.SS.15 


TRIDOOSO 


UO  4b  K=A,N  TRI00060 

AKSA-1  TR1D0070 


sum  contains  the  sum  of  the  souare  Elements 
of  a  column,  except  the  first  k-2  elements . 


SoM=A(A-1.K-2)»A(K-1,K-2)  thioooso 
UO  20  J=K,N  TR 100090 
SUM=SUM*A ( J » K-2  > • A ( J . K-2 )  TRI00100 


THE.  b  ARRAY  l0nTA1n$  THE  BETA  VALUES. 
(I.t.r  bETA=SUM«»(l/2) ) 


6<K-2l=obIGN(DSaRT(bUM) .-AlK-l.K-2) ) 


if  betaso  no  transformation  is  initiated. 

IF(B<Ko*>)  24,46.24  THID0120 


THE  COMPONENTS  OF  THE  COLUMN  vector  6  are 
STORED  IN  ThwRuSmONS  OF  THE  ANNIHILATED 
ELEMENTS  OF  a  (I.E.,  LOVER  HALF  OF  A) 


A(K-l.K-E)sUSeRTtu.S*CAdS<AO<-I,K«2)/a(K'2>  >«0.S> 


UENOM=-*.»AlA-l.R-2>4BlK-2>  TRI00140 
UO  30  |:inn  TRID0160 
AtI.A-2)=A(I,K-2)/UENOM  TRI00170 
bCALSU.  TRID0175 
UO  36  UbRK.N  TrIDOIGO 


THc.  BETAS  ARl  FORMED  ONE  BY  Of*. 

•HEN  R  OF  Th«.4  HAVE  BFEN  FORMED. 

V  ANL  9  contain  only  (N-R)  ELEMENTS. 
lMt  t  ARRAY  is  UED  TO  STORt  SUCCESSIVE 
P.b  *  MJ  0,S. 


JIUISO. 

iHb.tO.MTIOO  TO  Abo 
UO  349  LbAK.oU 
o<U)S4.(wt*A(U,LI«AtL,4'2) 


TRlr0190 

TR1C0191 

THID0192 

TRIP0I93 
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Rt.LY.GAUSsNE.lU> 

1 


SUB? 

SUBROUTINE  E16VAL(LP,E.A.B.«.F) 
1MPLIC1I  REAL.SIA-H.O-Z) 


10 

c 

11 

c 

Id 

c 

1  b 

c 

14 

c 

10 

UlMtNSIoN  E(| 

lu 

uiMLNSIu..  Fd 

17 

c 

Iti 

c 

iy 

c 

4.0 

Ai'irLALS  t  a  ( 1 ) 

*1 

LF=U. 

4 id 

Lw  I  1S.-.LP 

ANsDMAXIIAK.i 

*4 

1 

bi-.rUMAXi  IriM.I 

*b 

bUSAK+UMVljR 

do 

UC  b  IA..LP 

cl 

c 

do 

c 

*y 

c 

bo 

c 

Jl 

c 

-id 

c 

oJ 

A  ( I ) -A ( a  >  / bL 

34 

b l I 1 =L ( i ) / UL 

3b 

HI)=-1.U 

JO 

6 

«(!)=!. v 

J7 

UU  SO  N=i«LP 

ob 

c 

3y 

c 

40 

l 

41 

c 

4  d 

c 

4b 

c 

44 

c 

43 

L 

Ho 

c 

4  7 

d 

IF ( («(K)-E|N 

Ho 

lo 

A=(*(K)*L(K) 

ny 

c 

30 

L 

3l 

c 

3* 

c 

3b 

c 

3H 

Swsl.u 

3  j 

F<1>=AU>-X 

LP  IS  the  size  of  aepay  a. 

E  IS  A  VECToa  of  LP  ELEMFNTS  WHICH  will  HOLD 
THE  EIGENVALUES  in  descending  absolute  ORDER. 

A  is  A  VECTOR  OF  LP  ELEMENTS  GIVING  THE  DIAGONAL 
ELEMENTS  OF  THE  TRIDIAGONAL  MATRIX, 
b  IS  A  VECTOR  OF  LP  ELEMENTS.  THE  LAST  LP  -  1 
X  IS  A  VECTOR  OF  LP  ELEMENTS  USED  FOR  TEMPORARY 
STORAGE • 

F  AS  A  VECTOr  OF  LP  ELEMENTS  USED  FOR  TEMPORARY 
STORAGE . 


F  li.u  AbSOLUlo  dOUN^  FDR  THE  EIGENVALUES 


THIS  LOOP  FORCES  ThE  EIGENVALUES  TO  lIE  BETWEEN 
plus  ARE  MIIK.S  CI.E.  THE  E  AND  V  VECTupS  ARE 
HESRlCTIVElY  LU*  Ai.r  HIGH  ESTIMATES  ID  ALL 
Tml  rluENVALs.Cs, 


FIno  THE  K-Tm  EIGE.. VALUE.  ALSO  LOW  AiiD  HIGH 
ESTImaILS  FOr  THE  n-H-ST  to  LP-TH  EIGENVALUES 
ARt  IMPROVED.  THE  EIGENVALUES  A»E  Ful'M,  IN 
ascending  oi.uEr. 

THL  A-TH  EIOlNVaLUe  IS  CONSIDERED  FOuUD  IF  ThE 
THE  HIGH  ANu  LOW  Pl-ACES  AGREF  TO  SEVEN  UECIMAL 
PLACES. 


X  lb  A  GUESS  FOR  THE  X-TH  EIGENVALUE.  COMpuTt 
NUMBER  OF  EIGENVALUES  EQUAL  OR  EXCEEuINg  *  PV 
USING  STURM  SEQUENCE  (OPTEGA'S  METhOu). 


EvAL 


EVAU 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

evai 

FVAL 


EVAL 


EVAL 

EVAL 


EvAL 

rVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 


EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 


0.12 


37 

102 

Si=-1.0 

bo 

N=0 

b9 

60  TO  luG 

©u 

1UH 

Slsl.O 

ol 

N-l 

Ot 

103 

00  120  i=2. U 

©b 

IFCBIIM  106 

OH 

1UO 

lF(BlI-D)  li 

03 

1J7 

IK  (UABSiF  ( I—, 

uo 

C 

o  7 

c 

DO 

0 

u9 

c 

7J 

c 

7 i 

111 

F(I-l)=r ll-l 

It 

MI-2!:MI-t 

7b 

lie 

F(I)=(Ail)-A 

7h 

OU  10  ll3 

73 

lib 

F(II=(All)-X 

7o 

Go  Tu  lib 

77 

114 

F(J)=(A(i)-X 

7o 

1*3 

SisSl 

/9 

IF(F  (I niio.; 

00 

lib 

SIsLSIGnCSI.I 

©1 

lF(Sl«S.ni7 

oc 

117 

n=n«  1 

03 

ltu 

CUnTINUc. 

oH 

C 

u3 

c 

DO 

c 

©7 

c 

uo 

n=lP-N 

t>9 

IF  CN.LT. h)  O' 

90 

c 

91 

L 

9* 

c 

93 

c 

9h 

it 

uo  IS  OiA.N 

9j 

13 

.col=x 

9o 

to 

t.=N»i 

9  f 

0 

9o 

0 

99 

w 

luD 

c 

loi 

0 

iut 

lF(LP.Ll.N)  1 

103 

1/0  26  Jill.LP 

1U* 

0 

IwJ 

C 

10U 

0 

10/ 

C 

loo 

If  (X  •  LlJI) 

109 

to 

t ( J)=X 

Hu 

GO  TO  6 

ill 

bo 

CONI  IHOt 

1  At 

0 

lib 

c 

IF  1  Mt  PREVIOUS  T wo  TERMS  OF  THE  STUhv  SLOuF.nCE 
SEKt  VERT  SM.Li-i  ThFT  ARE  FORCED  TO  t)E  CLOSER 
TO  ONE  IN  MAGNITUDE  TO  AVOID  UNDERFLOW  PROBLEMS 


NO.  LlT  N  BL  TtiE  NU'TFR  OF 
E  Iotf. VALUES  sM.LLE.  Twjr:  t . 


x  otLOMES  a,,  upper  rouNc  for  the 
A-n.  To  r.-Th  eigenvalues. 


IF  all  THE  EIGENVALUES  ARE  SMALLER  TnAN 
TES1  .HE Then  Wl  HAVE  CONVERGED  TO  Thl 
K-Th  EIGENVALUl. 


if  a  IS  LARolR  than  PREVIOUS  LOWrR 
bOul.L.  INChlaSL  THl  LOwER  bound. 


RESloPt  INPol  AND  SCALE  EIGENVALUES. 


EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 


EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

EVAL 

FVAL 

EVAL 

EVAL 

EVAL 

FVAL 

FVAL 

EVAL 

EvaL 

Eval 

EVAL 

F.VAL 

FVAL 

Eval 

Eval 

eval 

eval 


eval 

fval 

EVAL 

EVAL 
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114 

116 

llo 

117 

llo 

119 

c 

60 
c _ _ 

UO  60  1=1 .LP 

A<I1=A(1(»BD 

BII)=BU)*BU 

»U>  =  (H1)+E<I>  >»6D*0.5 

EVAL 

EVAL 

EVAL 

Eval 

120 

c 

SOkl  EIGENVALUES  IN  ABSOLUTE  DESCENDING  ORDER 

EVAL 

122 

J=LP 

EVAL 

123 

K=l 

Eval 

124 

DO  SO  1=1 >Lk 

Fval 

l«:b 

1F(DABSI»(K) )-DAbS(»(J)))e3>63>65 

12o 

03 

E<l)=k(U> 

Eval 

U7 

J=J-1 

Eval 

12b 

GO  TO  80 

Eval 

129 

6b 

UIJSftOt) 

Eval 

loO 

k=k  +  1 

eval 

lii 

00 

CunT  11. Ul 

Eval 

102 

KLTOHN 

Eval 

lb3 

End 
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RllT.GAUSsnE.T  ( 1 ) 
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«♦ 
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b 

c 
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c 
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c 

a 

c 
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c 

10 

c 

11 

c 
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c 

1** 

c 

lb 

lu 

17 

lc 

IV 

dU 

c 

21 
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The  above  sample  output  listing  is  given  here  only  to  illustrate 
the  format  as  programmed  in  RELY  I  at  the  time  of  delivery.  The  values 
shown  are  of  no  significance  relative  to  the  content  of  the  report 
(though  they  are  those  of  a  sample  from  Example  IV,  Sec.  3.3),  nor  are 
they  expected  to  be  repeatable  exactly  with  implementation  of  RELY  I  at 
a  different  computer  tacility.  The  sample  is  offered  as  an  aid  to 
users,  showing  format  and  exemplary  behavior  in  a  given  random  case. 
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